A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
tensor.f90
Go to the documentation of this file.
1 
2 ! tensor.f90
3 !
4 !> Tensor utility module
5 !
6 !> @copyright
7 !> 2015-2018 Lesley De Cruz & Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !> coolist is a type and also means "coordinate list"
14 !
15 !---------------------------------------------------------------------------!
16 
17 
18 MODULE tensor
19  USE params, only: ndim
20  IMPLICIT NONE
21 
22  PRIVATE
23 
24  !> Coordinate list element type. Elementary elements of the sparse tensors.
25  TYPE :: coolist_elem
26  INTEGER :: j !< Index \f$j\f$ of the element
27  INTEGER :: k !< Index \f$k\f$ of the element
28  REAL(KIND=8) :: v !< Value of the element
29  END TYPE coolist_elem
30 
31  !> 4d coordinate list element type. Elementary elements of the 4d sparse tensors.
32  TYPE :: coolist_elem4
33  INTEGER :: j,k,l
34  REAL(KIND=8) :: v
35  END TYPE coolist_elem4
36 
37  !> Coordinate list. Type used to represent the sparse tensor.
38  TYPE, PUBLIC :: coolist
39  TYPE(coolist_elem), DIMENSION(:), ALLOCATABLE :: elems !< Lists of elements tensor::coolist_elem
40  INTEGER :: nelems = 0 !< Number of elements in the list.
41  END TYPE coolist
42 
43  !> 4d coordinate list. Type used to represent the rank-4 sparse tensor.
44  TYPE, PUBLIC :: coolist4
45  TYPE(coolist_elem4), DIMENSION(:), ALLOCATABLE :: elems
46  INTEGER :: nelems = 0
47  END TYPE coolist4
48 
49  !> Parameter to test the equality with zero.
50  REAL(KIND=8), PARAMETER :: real_eps = 2.2204460492503131e-16
51 
54 
57  PUBLIC :: print_tensor4
61  PUBLIC :: tensor_empty,tensor4_empty
64 
65 CONTAINS
66 
67  !> Routine to copy a coolist.
68  !> @param src Source coolist
69  !> @param dst Destination coolist
70  !> @remark The destination tensor have to be an empty tensor, i.e. with unallocated list of elements and nelems set to 0.
71  SUBROUTINE copy_coo(src,dst)
72  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: src
73  TYPE(coolist), DIMENSION(ndim), INTENT(OUT) :: dst
74  INTEGER :: i,j,AllocStat
75 
76  DO i=1,ndim
77  IF (dst(i)%nelems/=0) stop "*** copy_coo : Destination coolist not empty ! ***"
78  ALLOCATE(dst(i)%elems(src(i)%nelems), stat=allocstat)
79  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
80  DO j=1,src(i)%nelems
81  dst(i)%elems(j)%j=src(i)%elems(j)%j
82  dst(i)%elems(j)%k=src(i)%elems(j)%k
83  dst(i)%elems(j)%v=src(i)%elems(j)%v
84  ENDDO
85  dst(i)%nelems=src(i)%nelems
86  ENDDO
87  END SUBROUTINE copy_coo
88 
89  !> Routine to convert a matrix to a tensor.
90  !> @param src Source matrix
91  !> @param dst Destination tensor
92  !> @remark The destination tensor have to be an empty tensor, i.e. with unallocated list of elements and nelems set to 0.
93  SUBROUTINE mat_to_coo(src,dst)
94  REAL(KIND=8), DIMENSION(0:ndim,0:ndim), INTENT(IN) :: src
95  TYPE(coolist), DIMENSION(ndim), INTENT(OUT) :: dst
96  INTEGER :: i,j,n,AllocStat
97  DO i=1,ndim
98  n=0
99  DO j=1,ndim
100  IF (abs(src(i,j))>real_eps) n=n+1
101  ENDDO
102  IF (n/=0) THEN
103  IF (dst(i)%nelems/=0) stop "*** mat_to_coo : Destination coolist not empty ! ***"
104  ALLOCATE(dst(i)%elems(n), stat=allocstat)
105  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
106  n=0
107  DO j=1,ndim
108  IF (abs(src(i,j))>real_eps) THEN
109  n=n+1
110  dst(i)%elems(n)%j=j
111  dst(i)%elems(n)%k=0
112  dst(i)%elems(n)%v=src(i,j)
113  ENDIF
114  ENDDO
115  ENDIF
116  dst(i)%nelems=n
117  ENDDO
118  END SUBROUTINE mat_to_coo
119 
120  !> Sparse multiplication of a tensor with two vectors: \f${\displaystyle \sum_{j,k=0}^{ndim}} \mathcal{T}_{i,j,k} \, a_j \,b_k\f$.
121  !> @param coolist_ijk a coordinate list (sparse tensor) of which index
122  !> 2 and 3 will be contracted.
123  !> @param arr_j the vector to be contracted with index 2 of coolist_ijk
124  !> @param arr_k the vector to be contracted with index 3 of coolist_ijk
125  !> @param res vector (buffer) to store the result of the contraction
126  !> @remark Note that it is NOT safe to pass `arr_j`/`arr_k` as a result buffer,
127  !> as this operation does multiple passes.
128  SUBROUTINE sparse_mul3(coolist_ijk, arr_j, arr_k, res)
129  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ijk
130  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_j, arr_k
131  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
132  INTEGER :: i,j,k,n
133  res=0.d0
134  DO i=1,ndim
135  DO n=1,coolist_ijk(i)%nelems
136  j=coolist_ijk(i)%elems(n)%j
137  k=coolist_ijk(i)%elems(n)%k
138  res(i) = res(i) + coolist_ijk(i)%elems(n)%v * arr_j(j)*arr_k(k)
139  END DO
140  END DO
141  END SUBROUTINE sparse_mul3
142 
143  !> Sparse multiplication of two tensors to determine the Jacobian:
144  !> \f[J_{i,j} = {\displaystyle \sum_{k=0}^{ndim}} \left( \mathcal{T}_{i,j,k} + \mathcal{T}_{i,k,j} \right) \, a_k.\f]
145  !> It's implemented slightly differently: for every \f$\mathcal{T}_{i,j,k}\f$, we add to \f$J_{i,j}\f$ as follows:
146  !> \f[J_{i,j} = J_{i,j} + \mathcal{T}_{i,j,k} \, a_k \\ J_{i,k} = J_{i,k} + \mathcal{T}_{i,j,k} \, a_j\f]
147  !> This version return a coolist (sparse tensor).
148  !> @param coolist_ijk a coordinate list (sparse tensor) of which index
149  !> 2 or 3 will be contracted.
150  !> @param arr_j the vector to be contracted with index 2 and then index 3 of ffi_coo_ijk
151  !> @param jcoo_ij a coolist (sparse tensor) to store the result of the contraction
152  SUBROUTINE jsparse_mul(coolist_ijk, arr_j, jcoo_ij)
153  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ijk
154  TYPE(coolist), DIMENSION(ndim), INTENT(OUT):: jcoo_ij
155  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_j
156  REAL(KIND=8) :: v
157  INTEGER :: i,j,k,n,nj,AllocStat
158  DO i=1,ndim
159  IF (jcoo_ij(i)%nelems/=0) stop "*** jsparse_mul : Destination coolist not empty ! ***"
160  nj=2*coolist_ijk(i)%nelems
161  ALLOCATE(jcoo_ij(i)%elems(nj), stat=allocstat)
162  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
163  nj=0
164  DO n=1,coolist_ijk(i)%nelems
165  j=coolist_ijk(i)%elems(n)%j
166  k=coolist_ijk(i)%elems(n)%k
167  v=coolist_ijk(i)%elems(n)%v
168  IF (j /=0) THEN
169  nj=nj+1
170  jcoo_ij(i)%elems(nj)%j=j
171  jcoo_ij(i)%elems(nj)%k=0
172  jcoo_ij(i)%elems(nj)%v=v*arr_j(k)
173  END IF
174 
175  IF (k /=0) THEN
176  nj=nj+1
177  jcoo_ij(i)%elems(nj)%j=k
178  jcoo_ij(i)%elems(nj)%k=0
179  jcoo_ij(i)%elems(nj)%v=v*arr_j(j)
180  END IF
181  END DO
182  jcoo_ij(i)%nelems=nj
183  END DO
184  END SUBROUTINE jsparse_mul
185 
186  !> Sparse multiplication of two tensors to determine the Jacobian:
187  !> \f[J_{i,j} = {\displaystyle \sum_{k=0}^{ndim}} \left( \mathcal{T}_{i,j,k} + \mathcal{T}_{i,k,j} \right) \, a_k.\f]
188  !> It's implemented slightly differently: for every \f$\mathcal{T}_{i,j,k}\f$, we add to \f$J_{i,j}\f$ as follows:
189  !> \f[J_{i,j} = J_{i,j} + \mathcal{T}_{i,j,k} \, a_k \\ J_{i,k} = J_{i,k} + \mathcal{T}_{i,j,k} \, a_j\f]
190  !> This version return a matrix.
191  !> @param coolist_ijk a coordinate list (sparse tensor) of which index
192  !> 2 or 3 will be contracted.
193  !> @param arr_j the vector to be contracted with index 2 and then index 3 of ffi_coo_ijk
194  !> @param jcoo_ij a matrix to store the result of the contraction
195  SUBROUTINE jsparse_mul_mat(coolist_ijk, arr_j, jcoo_ij)
196  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ijk
197  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT):: jcoo_ij
198  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_j
199  REAL(KIND=8) :: v
200  INTEGER :: i,j,k,n
201  jcoo_ij=0.d0
202  DO i=1,ndim
203  DO n=1,coolist_ijk(i)%nelems
204  j=coolist_ijk(i)%elems(n)%j
205  k=coolist_ijk(i)%elems(n)%k
206  v=coolist_ijk(i)%elems(n)%v
207  IF (j /=0) jcoo_ij(i,j)=jcoo_ij(i,j)+v*arr_j(k)
208  IF (k /=0) jcoo_ij(i,k)=jcoo_ij(i,k)+v*arr_j(j)
209  END DO
210  END DO
211  END SUBROUTINE jsparse_mul_mat
212 
213  !> Sparse multiplication of a 2d sparse tensor with a vector: \f${\displaystyle \sum_{j=0}^{ndim}} \mathcal{T}_{i,j,k} \, a_j \f$.
214  !> @param coolist_ij a coordinate list (sparse tensor) of which index
215  !> 2 will be contracted.
216  !> @param arr_j the vector to be contracted with index 2 of coolist_ijk
217  !> @param res vector (buffer) to store the result of the contraction
218  !> @remark Note that it is NOT safe to pass `arr_j` as a result buffer,
219  !> as this operation does multiple passes.
220  SUBROUTINE sparse_mul2(coolist_ij, arr_j, res)
221  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ij
222  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_j
223  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
224  INTEGER :: i,j,n
225  res=0.d0
226  DO i=1,ndim
227  DO n=1,coolist_ij(i)%nelems
228  j=coolist_ij(i)%elems(n)%j
229  res(i) = res(i) + coolist_ij(i)%elems(n)%v * arr_j(j)
230  END DO
231  END DO
232  END SUBROUTINE sparse_mul2
233 
234  !> Routine to simplify a coolist (sparse tensor). For each index \f$i\f$, it upper triangularize the matrix
235  !> \f[\mathcal{T}_{i,j,k} \qquad 0 \leq j,k \leq ndim.\f]
236  !> @param tensor a coordinate list (sparse tensor) which will be simplified.
237  SUBROUTINE simplify(tensor)
238  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT):: tensor
239  INTEGER :: i,j,k
240  INTEGER :: li,lii,liii,n
241  DO i= 1,ndim
242  n=tensor(i)%nelems
243  DO li=n,2,-1
244  j=tensor(i)%elems(li)%j
245  k=tensor(i)%elems(li)%k
246  DO lii=li-1,1,-1
247  IF (((j==tensor(i)%elems(lii)%j).AND.(k==tensor(i)&
248  &%elems(lii)%k)).OR.((j==tensor(i)%elems(lii)%k).AND.(k==tensor(i)%elems(lii)%j))) THEN
249  ! Found another entry with the same i,j,k: merge both into
250  ! the one listed first (of those two).
251  tensor(i)%elems(lii)%v=tensor(i)%elems(lii)%v+tensor(i)%elems(li)%v
252  IF (j>k) THEN
253  tensor(i)%elems(lii)%j=tensor(i)%elems(li)%k
254  tensor(i)%elems(lii)%k=tensor(i)%elems(li)%j
255  ENDIF
256 
257  ! Shift the rest of the items one place down.
258  DO liii=li+1,n
259  tensor(i)%elems(liii-1)%j=tensor(i)%elems(liii)%j
260  tensor(i)%elems(liii-1)%k=tensor(i)%elems(liii)%k
261  tensor(i)%elems(liii-1)%v=tensor(i)%elems(liii)%v
262  END DO
263  tensor(i)%nelems=tensor(i)%nelems-1
264  ! Here we should stop because the li no longer points to the
265  ! original i,j,k element
266  EXIT
267  ENDIF
268  ENDDO
269  ENDDO
270  n=tensor(i)%nelems
271  DO li=1,n
272  ! Clear new "almost" zero entries and shift rest of the items one place down.
273  ! Make sure not to skip any entries while shifting!
274  DO WHILE (abs(tensor(i)%elems(li)%v) < real_eps)
275  DO liii=li+1,n
276  tensor(i)%elems(liii-1)%j=tensor(i)%elems(liii)%j
277  tensor(i)%elems(liii-1)%k=tensor(i)%elems(liii)%k
278  tensor(i)%elems(liii-1)%v=tensor(i)%elems(liii)%v
279  ENDDO
280  tensor(i)%nelems=tensor(i)%nelems-1
281  if (li > tensor(i)%nelems) THEN
282  EXIT
283  ENDIF
284  ENDDO
285  ENDDO
286 
287  n=tensor(i)%nelems
288  DO li=1,n
289  ! Upper triangularize
290  j=tensor(i)%elems(li)%j
291  k=tensor(i)%elems(li)%k
292  IF (j>k) THEN
293  tensor(i)%elems(li)%j=k
294  tensor(i)%elems(li)%k=j
295  ENDIF
296  ENDDO
297 
298 
299  ENDDO
300  END SUBROUTINE simplify
301 
302 
303  !> Subroutine to add element to a coolist.
304  !> @param t destination tensor
305  !> @param i tensor \f$i\f$ index
306  !> @param j tensor \f$j\f$ index
307  !> @param k tensor \f$k\f$ index
308  !> @param v value to add
309  SUBROUTINE add_elem(t,i,j,k,v)
310  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: t
311  INTEGER, INTENT(IN) :: i,j,k
312  REAL(KIND=8), INTENT(IN) :: v
313  INTEGER :: n
314  IF (abs(v) .ge. real_eps) THEN
315  n=(t(i)%nelems)+1
316  t(i)%elems(n)%j=j
317  t(i)%elems(n)%k=k
318  t(i)%elems(n)%v=v
319  t(i)%nelems=n
320  END IF
321  END SUBROUTINE add_elem
322 
323  !> Subroutine to add element to a coolist and check for overflow.
324  !> Once the t buffer tensor is full, add it to the destination buffer.
325  !> @param t temporary buffer tensor for the destination tensor
326  !> @param i tensor \f$i\f$ index
327  !> @param j tensor \f$j\f$ index
328  !> @param k tensor \f$k\f$ index
329  !> @param v value to add
330  !> @param dst destination tensor
331  SUBROUTINE add_check(t,i,j,k,v,dst)
332  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: t
333  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: dst
334  INTEGER, INTENT(IN) :: i,j,k
335  REAL(KIND=8), INTENT(IN) :: v
336  INTEGER :: n
337  CALL add_elem(t,i,j,k,v)
338  IF (t(i)%nelems==size(t(i)%elems)) THEN
339  CALL add_to_tensor(t,dst)
340  DO n=1,ndim
341  t(n)%nelems=0
342  ENDDO
343  ENDIF
344  END SUBROUTINE add_check
345 
346  !> Routine to add a rank-3 tensor to another one.
347  !> @param src Tensor to add
348  !> @param dst Destination tensor
349  SUBROUTINE add_to_tensor(src,dst)
350  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: src
351  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: dst
352  TYPE(coolist_elem), DIMENSION(:), ALLOCATABLE :: celems
353  INTEGER :: i,j,n,AllocStat
354 
355  DO i=1,ndim
356  IF (src(i)%nelems/=0) THEN
357  IF (dst(i)%nelems==0) THEN
358  IF (ALLOCATED(dst(i)%elems)) THEN
359  DEALLOCATE(dst(i)%elems, stat=allocstat)
360  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
361  ENDIF
362  ALLOCATE(dst(i)%elems(src(i)%nelems), stat=allocstat)
363  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
364  n=0
365  ELSE
366  n=dst(i)%nelems
367  ALLOCATE(celems(n), stat=allocstat)
368  DO j=1,n
369  celems(j)%j=dst(i)%elems(j)%j
370  celems(j)%k=dst(i)%elems(j)%k
371  celems(j)%v=dst(i)%elems(j)%v
372  ENDDO
373  DEALLOCATE(dst(i)%elems, stat=allocstat)
374  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
375  ALLOCATE(dst(i)%elems(src(i)%nelems+n), stat=allocstat)
376  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
377  DO j=1,n
378  dst(i)%elems(j)%j=celems(j)%j
379  dst(i)%elems(j)%k=celems(j)%k
380  dst(i)%elems(j)%v=celems(j)%v
381  ENDDO
382  DEALLOCATE(celems, stat=allocstat)
383  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
384  ENDIF
385  DO j=1,src(i)%nelems
386  dst(i)%elems(n+j)%j=src(i)%elems(j)%j
387  dst(i)%elems(n+j)%k=src(i)%elems(j)%k
388  dst(i)%elems(n+j)%v=src(i)%elems(j)%v
389  ENDDO
390  dst(i)%nelems=src(i)%nelems+n
391  ENDIF
392  ENDDO
393 
394  END SUBROUTINE add_to_tensor
395 
396  !> Routine to print a rank 3 tensor coolist.
397  !> @param t coolist to print
398  SUBROUTINE print_tensor(t,s)
399  USE util, only: str
400  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: t
401  CHARACTER, INTENT(IN), OPTIONAL :: s
402  CHARACTER :: r
403  INTEGER :: i,n,j,k
404  IF (PRESENT(s)) THEN
405  r=s
406  ELSE
407  r="t"
408  END IF
409  DO i=1,ndim
410  DO n=1,t(i)%nelems
411  j=t(i)%elems(n)%j
412  k=t(i)%elems(n)%k
413  IF( abs(t(i)%elems(n)%v) .GE. real_eps) THEN
414  write(*,"(A,ES12.5)") r//"["//trim(str(i))//"]["//trim(str(j)) &
415  &//"]["//trim(str(k))//"] = ",t(i)%elems(n)%v
416  END IF
417  END DO
418  END DO
419  END SUBROUTINE print_tensor
420 
421  !> Load a rank-4 tensor coolist from a file definition
422  !> @param s Destination filename
423  !> @param t The coolist to write
424  SUBROUTINE write_tensor_to_file(s,t)
425  CHARACTER (LEN=*), INTENT(IN) :: s
426  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: t
427  INTEGER :: i,j,k,n
428  OPEN(30,file=s)
429  DO i=1,ndim
430  WRITE(30,*) i,t(i)%nelems
431  DO n=1,t(i)%nelems
432  j=t(i)%elems(n)%j
433  k=t(i)%elems(n)%k
434  WRITE(30,*) i,j,k,t(i)%elems(n)%v
435  END DO
436  END DO
437  CLOSE(30)
438  END SUBROUTINE write_tensor_to_file
439 
440  !> Load a rank-4 tensor coolist from a file definition
441  !> @param s Filename of the tensor definition file
442  !> @param t The loaded coolist
443  !> @remark The destination tensor have to be an empty tensor, i.e. with unallocated list of elements and nelems set to 0.
444  SUBROUTINE load_tensor_from_file(s,t)
445  CHARACTER (LEN=*), INTENT(IN) :: s
446  TYPE(coolist), DIMENSION(ndim), INTENT(OUT) :: t
447  INTEGER :: i,ir,j,k,n,AllocStat
448  REAL(KIND=8) :: v
449  OPEN(30,file=s,status='old')
450  DO i=1,ndim
451  READ(30,*) ir,n
452  IF (n /= 0) THEN
453  ALLOCATE(t(i)%elems(n), stat=allocstat)
454  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
455  t(i)%nelems=n
456  ENDIF
457  DO n=1,t(i)%nelems
458  READ(30,*) ir,j,k,v
459  t(i)%elems(n)%j=j
460  t(i)%elems(n)%k=k
461  t(i)%elems(n)%v=v
462  ENDDO
463  END DO
464  CLOSE(30)
465  END SUBROUTINE load_tensor_from_file
466 
467 
468 
469  !> Routine to add a matrix to a rank-3 tensor.
470  !> @param i Add to tensor component i
471  !> @param src Matrix to add
472  !> @param dst Destination tensor
473  SUBROUTINE add_matc_to_tensor(i,src,dst)
474  INTEGER, INTENT(IN) :: i
475  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(IN) :: src
476  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: dst
477  TYPE(coolist_elem), DIMENSION(:), ALLOCATABLE :: celems
478  INTEGER :: j,k,r,n,nsrc,AllocStat
479 
480  nsrc=0
481  DO j=1,ndim
482  DO k=1,ndim
483  IF (abs(src(j,k))>real_eps) nsrc=nsrc+1
484  END DO
485  END DO
486 
487  IF (dst(i)%nelems==0) THEN
488  IF (ALLOCATED(dst(i)%elems)) THEN
489  DEALLOCATE(dst(i)%elems, stat=allocstat)
490  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
491  ENDIF
492  ALLOCATE(dst(i)%elems(nsrc), stat=allocstat)
493  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
494  n=0
495  ELSE
496  n=dst(i)%nelems
497  ALLOCATE(celems(n), stat=allocstat)
498  DO j=1,n
499  celems(j)%j=dst(i)%elems(j)%j
500  celems(j)%k=dst(i)%elems(j)%k
501  celems(j)%v=dst(i)%elems(j)%v
502  ENDDO
503  DEALLOCATE(dst(i)%elems, stat=allocstat)
504  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
505  ALLOCATE(dst(i)%elems(nsrc+n), stat=allocstat)
506  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
507  DO j=1,n
508  dst(i)%elems(j)%j=celems(j)%j
509  dst(i)%elems(j)%k=celems(j)%k
510  dst(i)%elems(j)%v=celems(j)%v
511  ENDDO
512  DEALLOCATE(celems, stat=allocstat)
513  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
514  ENDIF
515  r=0
516  DO j=1,ndim
517  DO k=1,ndim
518  IF (abs(src(j,k))>real_eps) THEN
519  r=r+1
520  dst(i)%elems(n+r)%j=j
521  dst(i)%elems(n+r)%k=k
522  dst(i)%elems(n+r)%v=src(j,k)
523  ENDIF
524  ENDDO
525  END DO
526  dst(i)%nelems=nsrc+n
527 
528  END SUBROUTINE add_matc_to_tensor
529 
530 
531  !> Routine to add a matrix to a rank-4 tensor.
532  !> @param i Add to tensor component i,j
533  !> @param j Add to tensor component i,j
534  !> @param src Matrix to add
535  !> @param dst Destination tensor
536  SUBROUTINE add_matc_to_tensor4(i,j,src,dst)
537  INTEGER, INTENT(IN) :: i,j
538  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(IN) :: src
539  TYPE(coolist4), DIMENSION(ndim), INTENT(INOUT) :: dst
540  TYPE(coolist_elem4), DIMENSION(:), ALLOCATABLE :: celems
541  INTEGER :: k,l,r,n,nsrc,AllocStat
542 
543  nsrc=0
544  DO k=1,ndim
545  DO l=1,ndim
546  IF (abs(src(k,l))>real_eps) nsrc=nsrc+1
547  END DO
548  END DO
549 
550  IF (dst(i)%nelems==0) THEN
551  IF (ALLOCATED(dst(i)%elems)) THEN
552  DEALLOCATE(dst(i)%elems, stat=allocstat)
553  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
554  ENDIF
555  ALLOCATE(dst(i)%elems(nsrc), stat=allocstat)
556  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
557  n=0
558  ELSE
559  n=dst(i)%nelems
560  ALLOCATE(celems(n), stat=allocstat)
561  DO k=1,n
562  celems(k)%j=dst(i)%elems(k)%j
563  celems(k)%k=dst(i)%elems(k)%k
564  celems(k)%l=dst(i)%elems(k)%l
565  celems(k)%v=dst(i)%elems(k)%v
566  ENDDO
567  DEALLOCATE(dst(i)%elems, stat=allocstat)
568  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
569  ALLOCATE(dst(i)%elems(nsrc+n), stat=allocstat)
570  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
571  DO k=1,n
572  dst(i)%elems(k)%j=celems(k)%j
573  dst(i)%elems(k)%k=celems(k)%k
574  dst(i)%elems(k)%l=celems(k)%l
575  dst(i)%elems(k)%v=celems(k)%v
576  ENDDO
577  DEALLOCATE(celems, stat=allocstat)
578  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
579  ENDIF
580  r=0
581  DO k=1,ndim
582  DO l=1,ndim
583  IF (abs(src(k,l))>real_eps) THEN
584  r=r+1
585  dst(i)%elems(n+r)%j=j
586  dst(i)%elems(n+r)%k=k
587  dst(i)%elems(n+r)%l=l
588  dst(i)%elems(n+r)%v=src(k,l)
589  ENDIF
590  ENDDO
591  END DO
592  dst(i)%nelems=nsrc+n
593 
594  END SUBROUTINE add_matc_to_tensor4
595 
596 
597  !> Routine to add a vector to a rank-3 tensor.
598  !> @param j,k Add to tensor component j and k
599  !> @param src Vector to add
600  !> @param dst Destination tensor
601  SUBROUTINE add_vec_jk_to_tensor(j,k,src,dst)
602  INTEGER, INTENT(IN) :: j,k
603  REAL(KIND=8), DIMENSION(ndim), INTENT(IN) :: src
604  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: dst
605  TYPE(coolist_elem), DIMENSION(:), ALLOCATABLE :: celems
606  INTEGER :: i,l,r,n,nsrc,AllocStat
607 
608  DO i=1,ndim
609  nsrc=0
610  IF (abs(src(i))>real_eps) nsrc=1
611  IF (dst(i)%nelems==0) THEN
612  IF (ALLOCATED(dst(i)%elems)) THEN
613  DEALLOCATE(dst(i)%elems, stat=allocstat)
614  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
615  ENDIF
616  ALLOCATE(dst(i)%elems(nsrc), stat=allocstat)
617  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
618  n=0
619  ELSE
620  n=dst(i)%nelems
621  ALLOCATE(celems(n), stat=allocstat)
622  DO l=1,n
623  celems(l)%j=dst(i)%elems(l)%j
624  celems(l)%k=dst(i)%elems(l)%k
625  celems(l)%v=dst(i)%elems(l)%v
626  ENDDO
627  DEALLOCATE(dst(i)%elems, stat=allocstat)
628  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
629  ALLOCATE(dst(i)%elems(nsrc+n), stat=allocstat)
630  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
631  DO l=1,n
632  dst(i)%elems(l)%j=celems(l)%j
633  dst(i)%elems(l)%k=celems(l)%k
634  dst(i)%elems(l)%v=celems(l)%v
635  ENDDO
636  DEALLOCATE(celems, stat=allocstat)
637  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
638  ENDIF
639  r=0
640  IF (abs(src(i))>real_eps) THEN
641  r=r+1
642  dst(i)%elems(n+r)%j=j
643  dst(i)%elems(n+r)%k=k
644  dst(i)%elems(n+r)%v=src(i)
645  ENDIF
646  dst(i)%nelems=nsrc+n
647  END DO
648 
649 
650  END SUBROUTINE add_vec_jk_to_tensor
651 
652  !> Routine to add a vector to a rank-4 tensor plus permutation.
653  !> @param i,k,l Add to tensor component i,k and l
654  !> @param src Vector to add
655  !> @param dst Destination tensor
656  SUBROUTINE add_vec_ikl_to_tensor4_perm(i,k,l,src,dst)
657  INTEGER, INTENT(IN) :: i,k,l
658  REAL(KIND=8), DIMENSION(ndim), INTENT(IN) :: src
659  TYPE(coolist4), DIMENSION(ndim), INTENT(INOUT) :: dst
660  TYPE(coolist_elem4), DIMENSION(:), ALLOCATABLE :: celems
661  INTEGER :: j,ne,r,n,nsrc,AllocStat
662 
663  nsrc=0
664  DO j=1,ndim
665  IF (abs(src(j))>real_eps) nsrc=nsrc+1
666  ENDDO
667  nsrc=nsrc*3
668  IF (dst(i)%nelems==0) THEN
669  IF (ALLOCATED(dst(i)%elems)) THEN
670  DEALLOCATE(dst(i)%elems, stat=allocstat)
671  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
672  ENDIF
673  ALLOCATE(dst(i)%elems(nsrc), stat=allocstat)
674  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
675  n=0
676  ELSE
677  n=dst(i)%nelems
678  ALLOCATE(celems(n), stat=allocstat)
679  DO ne=1,n
680  celems(ne)%j=dst(i)%elems(ne)%j
681  celems(ne)%k=dst(i)%elems(ne)%k
682  celems(ne)%l=dst(i)%elems(ne)%l
683  celems(ne)%v=dst(i)%elems(ne)%v
684  ENDDO
685  DEALLOCATE(dst(i)%elems, stat=allocstat)
686  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
687  ALLOCATE(dst(i)%elems(nsrc+n), stat=allocstat)
688  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
689  DO ne=1,n
690  dst(i)%elems(ne)%j=celems(ne)%j
691  dst(i)%elems(ne)%k=celems(ne)%k
692  dst(i)%elems(ne)%l=celems(ne)%l
693  dst(i)%elems(ne)%v=celems(ne)%v
694  ENDDO
695  DEALLOCATE(celems, stat=allocstat)
696  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
697  ENDIF
698  r=0
699  DO j=1,ndim
700  IF (abs(src(j))>real_eps) THEN
701  r=r+1
702  dst(i)%elems(n+r)%j=j
703  dst(i)%elems(n+r)%k=k
704  dst(i)%elems(n+r)%l=l
705  dst(i)%elems(n+r)%v=src(j)
706  r=r+1
707  dst(i)%elems(n+r)%j=k
708  dst(i)%elems(n+r)%k=l
709  dst(i)%elems(n+r)%l=j
710  dst(i)%elems(n+r)%v=src(j)
711  r=r+1
712  dst(i)%elems(n+r)%j=l
713  dst(i)%elems(n+r)%k=j
714  dst(i)%elems(n+r)%l=k
715  dst(i)%elems(n+r)%v=src(j)
716  ENDIF
717  ENDDO
718  dst(i)%nelems=nsrc+n
719  END SUBROUTINE add_vec_ikl_to_tensor4_perm
720 
721  !> Routine to add a vector to a rank-4 tensor.
722  !> @param i,k,l Add to tensor component i,k and l
723  !> @param src Vector to add
724  !> @param dst Destination tensor
725  SUBROUTINE add_vec_ikl_to_tensor4(i,k,l,src,dst)
726  INTEGER, INTENT(IN) :: i,k,l
727  REAL(KIND=8), DIMENSION(ndim), INTENT(IN) :: src
728  TYPE(coolist4), DIMENSION(ndim), INTENT(INOUT) :: dst
729  TYPE(coolist_elem4), DIMENSION(:), ALLOCATABLE :: celems
730  INTEGER :: j,ne,r,n,nsrc,AllocStat
731 
732  nsrc=0
733  DO j=1,ndim
734  IF (abs(src(j))>real_eps) nsrc=nsrc+1
735  ENDDO
736 
737  IF (dst(i)%nelems==0) THEN
738  IF (ALLOCATED(dst(i)%elems)) THEN
739  DEALLOCATE(dst(i)%elems, stat=allocstat)
740  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
741  ENDIF
742  ALLOCATE(dst(i)%elems(nsrc), stat=allocstat)
743  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
744  n=0
745  ELSE
746  n=dst(i)%nelems
747  ALLOCATE(celems(n), stat=allocstat)
748  DO ne=1,n
749  celems(ne)%j=dst(i)%elems(ne)%j
750  celems(ne)%k=dst(i)%elems(ne)%k
751  celems(ne)%l=dst(i)%elems(ne)%l
752  celems(ne)%v=dst(i)%elems(ne)%v
753  ENDDO
754  DEALLOCATE(dst(i)%elems, stat=allocstat)
755  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
756  ALLOCATE(dst(i)%elems(nsrc+n), stat=allocstat)
757  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
758  DO ne=1,n
759  dst(i)%elems(ne)%j=celems(ne)%j
760  dst(i)%elems(ne)%k=celems(ne)%k
761  dst(i)%elems(ne)%l=celems(ne)%l
762  dst(i)%elems(ne)%v=celems(ne)%v
763  ENDDO
764  DEALLOCATE(celems, stat=allocstat)
765  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
766  ENDIF
767  r=0
768  DO j=1,ndim
769  IF (abs(src(j))>real_eps) THEN
770  r=r+1
771  dst(i)%elems(n+r)%j=j
772  dst(i)%elems(n+r)%k=k
773  dst(i)%elems(n+r)%l=l
774  dst(i)%elems(n+r)%v=src(j)
775  ENDIF
776  ENDDO
777  dst(i)%nelems=nsrc+n
778  END SUBROUTINE add_vec_ikl_to_tensor4
779 
780  !> Routine to add a vector to a rank-4 tensor.
781  !> @param i,j,k Add to tensor component i,j and k
782  !> @param src Vector to add
783  !> @param dst Destination tensor
784  SUBROUTINE add_vec_ijk_to_tensor4(i,j,k,src,dst)
785  INTEGER, INTENT(IN) :: i,j,k
786  REAL(KIND=8), DIMENSION(ndim), INTENT(IN) :: src
787  TYPE(coolist4), DIMENSION(ndim), INTENT(INOUT) :: dst
788  TYPE(coolist_elem4), DIMENSION(:), ALLOCATABLE :: celems
789  INTEGER :: l,ne,r,n,nsrc,AllocStat
790 
791  nsrc=0
792  DO l=1,ndim
793  IF (abs(src(l))>real_eps) nsrc=nsrc+1
794  ENDDO
795 
796  IF (dst(i)%nelems==0) THEN
797  IF (ALLOCATED(dst(i)%elems)) THEN
798  DEALLOCATE(dst(i)%elems, stat=allocstat)
799  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
800  ENDIF
801  ALLOCATE(dst(i)%elems(nsrc), stat=allocstat)
802  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
803  n=0
804  ELSE
805  n=dst(i)%nelems
806  ALLOCATE(celems(n), stat=allocstat)
807  DO ne=1,n
808  celems(ne)%j=dst(i)%elems(ne)%j
809  celems(ne)%k=dst(i)%elems(ne)%k
810  celems(ne)%l=dst(i)%elems(ne)%l
811  celems(ne)%v=dst(i)%elems(ne)%v
812  ENDDO
813  DEALLOCATE(dst(i)%elems, stat=allocstat)
814  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
815  ALLOCATE(dst(i)%elems(nsrc+n), stat=allocstat)
816  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
817  DO ne=1,n
818  dst(i)%elems(ne)%j=celems(ne)%j
819  dst(i)%elems(ne)%k=celems(ne)%k
820  dst(i)%elems(ne)%l=celems(ne)%l
821  dst(i)%elems(ne)%v=celems(ne)%v
822  ENDDO
823  DEALLOCATE(celems, stat=allocstat)
824  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
825  ENDIF
826  r=0
827  DO l=1,ndim
828  IF (abs(src(l))>real_eps) THEN
829  r=r+1
830  dst(i)%elems(n+r)%j=j
831  dst(i)%elems(n+r)%k=k
832  dst(i)%elems(n+r)%l=l
833  dst(i)%elems(n+r)%v=src(l)
834  ENDIF
835  ENDDO
836  dst(i)%nelems=nsrc+n
837  END SUBROUTINE add_vec_ijk_to_tensor4
838 
839 
840 
841 
842  !> Routine to convert a rank-3 tensor from matrix to coolist representation.
843  !> @param src Source matrix
844  !> @param dst Destination coolist
845  !> @remark The destination coolist have to be an empty one, i.e. with unallocated list of elements and nelems set to 0.
846  SUBROUTINE tensor_to_coo(src,dst)
847  REAL(KIND=8), DIMENSION(ndim,0:ndim,0:ndim), INTENT(IN) :: src
848  TYPE(coolist), DIMENSION(ndim), INTENT(OUT) :: dst
849  INTEGER :: i,j,k,n,AllocStat
850 
851  DO i=1,ndim
852  n=0
853  DO j=0,ndim
854  DO k=0,ndim
855  IF (abs(src(i,j,k))>real_eps) n=n+1
856  ENDDO
857  ENDDO
858  IF (n/=0) THEN
859  IF (dst(i)%nelems/=0) stop "*** tensor_to_coo : Destination coolist not empty ! ***"
860  ALLOCATE(dst(i)%elems(n), stat=allocstat)
861  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
862  n=0
863  DO j=0,ndim
864  DO k=0,ndim
865  IF (abs(src(i,j,k))>real_eps) THEN
866  n=n+1
867  dst(i)%elems(n)%j=j
868  dst(i)%elems(n)%k=k
869  dst(i)%elems(n)%v=src(i,j,k)
870  ENDIF
871  ENDDO
872  ENDDO
873  ENDIF
874  dst(i)%nelems=n
875  ENDDO
876  END SUBROUTINE tensor_to_coo
877 
878  !> Routine to convert a rank-4 tensor from matrix to coolist representation.
879  !> @param src Source matrix
880  !> @param dst Destination coolist
881  !> @remark The destination coolist have to be an empty one, i.e. with unallocated list of elements and nelems set to 0.
882  SUBROUTINE tensor4_to_coo4(src,dst)
883  REAL(KIND=8), DIMENSION(ndim,0:ndim,0:ndim,0:ndim), INTENT(IN) :: src
884  TYPE(coolist4), DIMENSION(ndim), INTENT(OUT) :: dst
885  INTEGER :: i,j,k,l,n,AllocStat
886 
887  DO i=1,ndim
888  n=0
889  DO j=0,ndim
890  DO k=0,ndim
891  DO l=0,ndim
892  IF (abs(src(i,j,k,l))>real_eps) n=n+1
893  ENDDO
894  ENDDO
895  ENDDO
896  IF (n/=0) THEN
897  IF (dst(i)%nelems/=0) stop "*** tensor_to_coo : Destination coolist not empty ! ***"
898  ALLOCATE(dst(i)%elems(n), stat=allocstat)
899  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
900  n=0
901  DO j=0,ndim
902  DO k=0,ndim
903  DO l=0,ndim
904  IF (abs(src(i,j,k,l))>real_eps) THEN
905  n=n+1
906  dst(i)%elems(n)%j=j
907  dst(i)%elems(n)%k=k
908  dst(i)%elems(n)%l=l
909  dst(i)%elems(n)%v=src(i,j,k,l)
910  ENDIF
911  ENDDO
912  ENDDO
913  ENDDO
914  ENDIF
915  dst(i)%nelems=n
916  ENDDO
917  END SUBROUTINE tensor4_to_coo4
918 
919  !> Routine to print a rank-4 tensor coolist.
920  !> @param t coolist to print
921  SUBROUTINE print_tensor4(t)
922  USE util, only: str
923  TYPE(coolist4), DIMENSION(ndim), INTENT(IN) :: t
924  INTEGER :: i,n,j,k,l
925  DO i=1,ndim
926  DO n=1,t(i)%nelems
927  j=t(i)%elems(n)%j
928  k=t(i)%elems(n)%k
929  l=t(i)%elems(n)%l
930  IF( abs(t(i)%elems(n)%v) .GE. real_eps) THEN
931  write(*,"(A,ES12.5)") "tensor["//trim(str(i))//"]["//trim(str(j)) &
932  &//"]["//trim(str(k))//"]["//trim(str(l))//"] = ",t(i)%elems(n)%v
933  END IF
934  END DO
935  END DO
936  END SUBROUTINE print_tensor4
937 
938 
939  !> Sparse multiplication of a rank-3 tensor coolist with a vector:
940  !> \f${\displaystyle \sum_{k=0}^{ndim}} \mathcal{T}_{i,j,k} \, b_k\f$.
941  !> Its output is a matrix.
942  !> @param coolist_ijk a coolist (sparse tensor) of which index k will be contracted.
943  !> @param arr_k the vector to be contracted with index k of coolist_ijk
944  !> @param res matrix (buffer) to store the result of the contraction
945  !> @remark Note that it is NOT safe to pass `arr_k` as a result buffer,
946  !> as this operation does multiple passes.
947  SUBROUTINE sparse_mul3_mat(coolist_ijk, arr_k, res)
948  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ijk
949  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_k
950  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: res
951  INTEGER :: i,j,k,n
952  res=0.d0
953  DO i=1,ndim
954  DO n=1,coolist_ijk(i)%nelems
955  j=coolist_ijk(i)%elems(n)%j
956  IF (j /= 0) THEN
957  k=coolist_ijk(i)%elems(n)%k
958  res(i,j) = res(i,j) + coolist_ijk(i)%elems(n)%v * arr_k(k)
959  ENDIF
960  END DO
961  END DO
962  END SUBROUTINE sparse_mul3_mat
963 
964 
965  !> Sparse multiplication of a rank-4 tensor coolist with three vectors: \f${\displaystyle \sum_{j,k,l=0}^{ndim}} \mathcal{T}_{i,j,k,l} \, a_j \,b_k \, c_l \f$.
966  !> @param coolist_ijkl a coolist (sparse tensor) of which index j, k and l will be contracted.
967  !> @param arr_j the vector to be contracted with index j of coolist_ijkl
968  !> @param arr_k the vector to be contracted with index k of coolist_ijkl
969  !> @param arr_l the vector to be contracted with index l of coolist_ijkl
970  !> @param res vector (buffer) to store the result of the contraction
971  !> @remark Note that it is NOT safe to pass `arr_j`/`arr_k`/`arr_l` as a result buffer,
972  !> as this operation does multiple passes.
973  SUBROUTINE sparse_mul4(coolist_ijkl, arr_j, arr_k, arr_l, res)
974  TYPE(coolist4), DIMENSION(ndim), INTENT(IN):: coolist_ijkl
975  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_j, arr_k, arr_l
976  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
977  INTEGER :: i,j,k,n,l
978  res=0.d0
979  DO i=1,ndim
980  DO n=1,coolist_ijkl(i)%nelems
981  j=coolist_ijkl(i)%elems(n)%j
982  k=coolist_ijkl(i)%elems(n)%k
983  l=coolist_ijkl(i)%elems(n)%l
984  res(i) = res(i) + coolist_ijkl(i)%elems(n)%v * arr_j(j)*arr_k(k)*arr_l(l)
985  END DO
986  END DO
987  END SUBROUTINE sparse_mul4
988 
989  !> Sparse multiplication of a tensor with two vectors: \f${\displaystyle \sum_{k,l=0}^{ndim}} \mathcal{T}_{i,j,k,l} \,b_k \, c_l \f$.
990  !> @param coolist_ijkl a coordinate list (sparse tensor) of which index
991  !> 3 and 4 will be contracted.
992  !> @param arr_k the vector to be contracted with index 3 of coolist_ijkl
993  !> @param arr_l the vector to be contracted with index 4 of coolist_ijkl
994  !> @param res matrix (buffer) to store the result of the contraction
995  !> @remark Note that it is NOT safe to pass `arr_k`/`arr_l` as a result buffer,
996  !> as this operation does multiple passes.
997  SUBROUTINE sparse_mul4_mat(coolist_ijkl, arr_k, arr_l, res)
998  TYPE(coolist4), DIMENSION(ndim), INTENT(IN):: coolist_ijkl
999  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_k, arr_l
1000  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: res
1001  INTEGER :: i,j,k,n,l
1002  res=0.d0
1003  DO i=1,ndim
1004  DO n=1,coolist_ijkl(i)%nelems
1005  j=coolist_ijkl(i)%elems(n)%j
1006  IF (j /= 0) THEN
1007  k=coolist_ijkl(i)%elems(n)%k
1008  l=coolist_ijkl(i)%elems(n)%l
1009  res(i,j) = res(i,j) + coolist_ijkl(i)%elems(n)%v * arr_k(k) * arr_l(l)
1010  ENDIF
1011  END DO
1012  END DO
1013  END SUBROUTINE sparse_mul4_mat
1014 
1015 
1016  !> Sparse multiplication of a 3d sparse tensor with a vectors: \f${\displaystyle \sum_{j=0}^{ndim}} \mathcal{T}_{i,j,k} \, a_j \f$.
1017  !> @param coolist_ijk a coordinate list (sparse tensor) of which index
1018  !> 2 will be contracted.
1019  !> @param arr_j the vector to be contracted with index 2 of coolist_ijk
1020  !> @param res vector (buffer) to store the result of the contraction
1021  !> @remark Note that it is NOT safe to pass `arr_j` as a result buffer,
1022  !> as this operation does multiple passes.
1023  SUBROUTINE sparse_mul2_j(coolist_ijk, arr_j, res)
1024  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ijk
1025  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_j
1026  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
1027  INTEGER :: i,j,n
1028  res=0.d0
1029  DO i=1,ndim
1030  DO n=1,coolist_ijk(i)%nelems
1031  j=coolist_ijk(i)%elems(n)%j
1032  res(i) = res(i) + coolist_ijk(i)%elems(n)%v * arr_j(j)
1033  END DO
1034  END DO
1035  END SUBROUTINE sparse_mul2_j
1036 
1037  !> Sparse multiplication of a rank-3 sparse tensor coolist with a vector: \f${\displaystyle \sum_{k=0}^{ndim}} \mathcal{T}_{i,j,k} \, a_k \f$.
1038  !> @param coolist_ijk a coordinate list (sparse tensor) of which index
1039  !> k will be contracted.
1040  !> @param arr_k the vector to be contracted with index k of coolist_ijk
1041  !> @param res vector (buffer) to store the result of the contraction
1042  !> @remark Note that it is NOT safe to pass `arr_k` as a result buffer,
1043  !> as this operation does multiple passes.
1044  SUBROUTINE sparse_mul2_k(coolist_ijk, arr_k, res)
1045  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ijk
1046  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: arr_k
1047  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
1048  INTEGER :: i,k,n
1049  res=0.d0
1050  DO i=1,ndim
1051  DO n=1,coolist_ijk(i)%nelems
1052  k=coolist_ijk(i)%elems(n)%k
1053  res(i) = res(i) + coolist_ijk(i)%elems(n)%v * arr_k(k)
1054  END DO
1055  END DO
1056  END SUBROUTINE sparse_mul2_k
1057 
1058 
1059  !> Routine to convert a rank-3 tensor coolist component into a matrix with i and k indices.
1060  !> @param src Source tensor
1061  !> @param dst Destination matrix
1062  SUBROUTINE coo_to_mat_ik(src,dst)
1063  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: src
1064  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: dst
1065  INTEGER :: i,n
1066 
1067  dst=0.d0
1068  DO i=1,ndim
1069  DO n=1,src(i)%nelems
1070  dst(i,src(i)%elems(n)%k)=src(i)%elems(n)%v
1071  ENDDO
1072  ENDDO
1073  END SUBROUTINE coo_to_mat_ik
1074 
1075  !> Routine to convert a rank-3 tensor coolist component into a matrix with i and j indices.
1076  !> @param src Source tensor
1077  !> @param dst Destination matrix
1078  SUBROUTINE coo_to_mat_ij(src,dst)
1079  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: src
1080  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: dst
1081  INTEGER :: i,n
1082 
1083  dst=0.d0
1084  DO i=1,ndim
1085  DO n=1,src(i)%nelems
1086  dst(i,src(i)%elems(n)%j)=src(i)%elems(n)%v
1087  ENDDO
1088  ENDDO
1089  END SUBROUTINE coo_to_mat_ij
1090 
1091  ! SUBROUTINE tensor_perm_ij(t)
1092  ! TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: t
1093  ! INTEGER :: i,j,k,n
1094 
1095  ! DO i=1,ndim
1096  ! DO n=1,t(i)%nelems
1097  ! j=t(i)%elems(n)%j
1098  ! k=t(i)%elems(n)%k
1099 
1100  ! t(i)%elems(n)%v
1101  ! ENDDO
1102  ! ENDDO
1103  ! END SUBROUTINE tensor_perm_ij
1104 
1105 !!! not so cool
1106 
1107  !> Routine to convert a rank-3 tensor coolist component into a matrix.
1108  !> @param i Component to convert
1109  !> @param src Source tensor
1110  !> @param dst Destination matrix
1111  SUBROUTINE coo_to_mat_i(i,src,dst)
1112  INTEGER, INTENT(IN) :: i
1113  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: src
1114  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: dst
1115  INTEGER :: n
1116 
1117  dst=0.d0
1118  DO n=1,src(i)%nelems
1119  dst(src(i)%elems(n)%j,src(i)%elems(n)%k)=src(i)%elems(n)%v
1120  ENDDO
1121  END SUBROUTINE coo_to_mat_i
1122 
1123  !> Routine to convert a rank-3 tensor coolist component into a vector.
1124  !> @param j Component j,k to convert
1125  !> @param k Component j,k to convert
1126  !> @param src Source tensor
1127  !> @param dst Destination vector
1128  SUBROUTINE coo_to_vec_jk(j,k,src,dst)
1129  INTEGER, INTENT(IN) :: j,k
1130  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: src
1131  REAL(KIND=8), DIMENSION(ndim), INTENT(OUT) :: dst
1132  INTEGER :: i,n
1133 
1134  dst=0.d0
1135  DO i=1,ndim
1136  DO n=1,src(i)%nelems
1137  IF ((src(i)%elems(n)%j==j).and.(src(i)%elems(n)%k==k)) dst(i)=src(i)%elems(n)%v
1138  END DO
1139  ENDDO
1140  END SUBROUTINE coo_to_vec_jk
1141 
1142 
1143  !> Routine to convert a rank-3 tensor coolist component into a matrix.
1144  !> @param j Component to convert
1145  !> @param src Source tensor
1146  !> @param dst Destination matrix
1147  SUBROUTINE coo_to_mat_j(j,src,dst)
1148  INTEGER, INTENT(IN) :: j
1149  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: src
1150  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: dst
1151  INTEGER :: i,n
1152 
1153  dst=0.d0
1154  DO i=1,ndim
1155  DO n=1,src(i)%nelems
1156  IF (src(i)%elems(n)%j==j) dst(i,src(i)%elems(n)%k)=src(i)%elems(n)%v
1157  ENDDO
1158  END DO
1159  END SUBROUTINE coo_to_mat_j
1160 
1161 
1162  !> Sparse multiplication of a rank-4 tensor coolist with a matrix : \f${\displaystyle \sum_{j,l=0}^{ndim}} \mathcal{T}_{i,j,k,l} \, m_{j,l} \f$.
1163  !> @param coolist_ijkl a coolist (sparse tensor) of which index j and l will be contracted.
1164  !> @param mat_jl the matrix to be contracted with indices j and l of coolist_ijkl
1165  !> @param res matrix (buffer) to store the result of the contraction
1166  !> @remark Note that it is NOT safe to pass `mat_jl` as a result buffer,
1167  !> as this operation does multiple passes.
1168  SUBROUTINE sparse_mul4_with_mat_jl(coolist_ijkl,mat_jl,res)
1169  TYPE(coolist4), DIMENSION(ndim), INTENT(IN):: coolist_ijkl
1170  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(IN) :: mat_jl
1171  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: res
1172  INTEGER i,j,k,l,n
1173 
1174  res=0.d0
1175  DO i=1,ndim
1176  DO n=1,coolist_ijkl(i)%nelems
1177  j=coolist_ijkl(i)%elems(n)%j
1178  k=coolist_ijkl(i)%elems(n)%k
1179  l=coolist_ijkl(i)%elems(n)%l
1180 
1181  res(i,k) = res(i,k) + coolist_ijkl(i)%elems(n)%v * mat_jl(j,l)
1182  ENDDO
1183  END DO
1184 
1185  END SUBROUTINE sparse_mul4_with_mat_jl
1186 
1187  !> Sparse multiplication of a rank-4 tensor coolist with a matrix : \f${\displaystyle \sum_{j,l=0}^{ndim}} \mathcal{T}_{i,j,k,l} \, m_{k,l} \f$.
1188  !> @param coolist_ijkl a coolist (sparse tensor) of which index k and l will be contracted.
1189  !> @param mat_kl the matrix to be contracted with indices k and l of coolist_ijkl
1190  !> @param res matrix (buffer) to store the result of the contraction
1191  !> @remark Note that it is NOT safe to pass `mat_kl` as a result buffer,
1192  !> as this operation does multiple passes.
1193  SUBROUTINE sparse_mul4_with_mat_kl(coolist_ijkl,mat_kl,res)
1194  TYPE(coolist4), DIMENSION(ndim), INTENT(IN):: coolist_ijkl
1195  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(IN) :: mat_kl
1196  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(OUT) :: res
1197  INTEGER i,j,k,l,n
1198 
1199  res=0.d0
1200  DO i=1,ndim
1201  DO n=1,coolist_ijkl(i)%nelems
1202  j=coolist_ijkl(i)%elems(n)%j
1203  k=coolist_ijkl(i)%elems(n)%k
1204  l=coolist_ijkl(i)%elems(n)%l
1205 
1206  res(i,j) = res(i,j) + coolist_ijkl(i)%elems(n)%v * mat_kl(k,l)
1207  ENDDO
1208  END DO
1209 
1210  END SUBROUTINE sparse_mul4_with_mat_kl
1211 
1212  !> Sparse multiplication of a rank-3 tensor coolist with a matrix: \f${\displaystyle \sum_{j,k=0}^{ndim}} \mathcal{T}_{i,j,k} \, m_{j,k}\f$.
1213  !> @param coolist_ijk a coolist (sparse tensor) of which index
1214  !> j and k will be contracted.
1215  !> @param mat_jk the matrix to be contracted with index j and k of coolist_ijk
1216  !> @param res vector (buffer) to store the result of the contraction
1217  !> @remark Note that it is NOT safe to pass `mat_jk` as a result buffer,
1218  !> as this operation does multiple passes.
1219  SUBROUTINE sparse_mul3_with_mat(coolist_ijk,mat_jk,res)
1220  TYPE(coolist), DIMENSION(ndim), INTENT(IN):: coolist_ijk
1221  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(IN) :: mat_jk
1222  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
1223  INTEGER i,j,k,n
1224 
1225  res=0.d0
1226  DO i=1,ndim
1227  DO n=1,coolist_ijk(i)%nelems
1228  j=coolist_ijk(i)%elems(n)%j
1229  k=coolist_ijk(i)%elems(n)%k
1230 
1231  res(i) = res(i) + coolist_ijk(i)%elems(n)%v * mat_jk(j,k)
1232  ENDDO
1233  END DO
1234 
1235  END SUBROUTINE sparse_mul3_with_mat
1236 
1237 
1238  !> Routine to convert a matrix to a rank-3 tensor.
1239  !> @param src Source matrix
1240  !> @param dst Destination tensor
1241  !> @remark The destination tensor have to be an empty tensor, i.e. with unallocated list of elements and nelems set to 0.
1242  !> @remark The j component will be set to 0.
1243  SUBROUTINE matc_to_coo(src,dst)
1244  REAL(KIND=8), DIMENSION(ndim,ndim), INTENT(IN) :: src
1245  TYPE(coolist), DIMENSION(ndim), INTENT(OUT) :: dst
1246  INTEGER :: i,j,n,AllocStat
1247  DO i=1,ndim
1248  n=0
1249  DO j=1,ndim
1250  IF (abs(src(i,j))>real_eps) n=n+1
1251  ENDDO
1252  IF (n/=0) THEN
1253  IF (dst(i)%nelems/=0) stop "*** mat_to_coo : Destination coolist not empty ! ***"
1254  ALLOCATE(dst(i)%elems(n), stat=allocstat)
1255  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
1256  n=0
1257  DO j=1,ndim
1258  IF (abs(src(i,j))>real_eps) THEN
1259  n=n+1
1260  dst(i)%elems(n)%j=0
1261  dst(i)%elems(n)%k=j
1262  dst(i)%elems(n)%v=src(i,j)
1263  ENDIF
1264  ENDDO
1265  ENDIF
1266  dst(i)%nelems=n
1267  ENDDO
1268  END SUBROUTINE matc_to_coo
1269 
1270  !> Routine to multiply a rank-3 tensor by a scalar
1271  !> @param s The scalar
1272  !> @param t The tensor
1273  SUBROUTINE scal_mul_coo(s,t)
1274  REAL(KIND=8), INTENT(IN) :: s
1275  TYPE(coolist), DIMENSION(ndim), INTENT(INOUT) :: t
1276  INTEGER :: i,li,n
1277  DO i=1,ndim
1278  n=t(i)%nelems
1279  DO li=1,n
1280  t(i)%elems(li)%v=s*t(i)%elems(li)%v
1281  ENDDO
1282  ENDDO
1283  END SUBROUTINE scal_mul_coo
1284 
1285  !> Test if a rank-3 tensor coolist is empty
1286  !> @param t rank-3 tensor coolist to be tested
1287  FUNCTION tensor_empty(t)
1288  TYPE(coolist), DIMENSION(ndim), INTENT(IN) :: t
1289  LOGICAL :: tensor_empty
1290  INTEGER :: i
1291  tensor_empty=.true.
1292  DO i=1,ndim
1293  IF (t(i)%nelems /= 0) THEN
1294  tensor_empty=.false.
1295  RETURN
1296  ENDIF
1297  END DO
1298  RETURN
1299  END FUNCTION tensor_empty
1300 
1301  !> Test if a rank-4 tensor coolist is empty
1302  !> @param t rank-4 tensor coolist to be tested
1303  FUNCTION tensor4_empty(t)
1304  TYPE(coolist4), DIMENSION(ndim), INTENT(IN) :: t
1305  LOGICAL :: tensor4_empty
1306  INTEGER :: i
1307  tensor4_empty=.true.
1308  DO i=1,ndim
1309  IF (t(i)%nelems /= 0) THEN
1310  tensor4_empty=.false.
1311  RETURN
1312  ENDIF
1313  END DO
1314  RETURN
1315  END FUNCTION tensor4_empty
1316 
1317  !> Load a rank-4 tensor coolist from a file definition
1318  !> @param s Filename of the tensor definition file
1319  !> @param t The loaded coolist
1320  !> @remark The destination tensor have to be an empty tensor, i.e. with unallocated list of elements and nelems set to 0.
1321  SUBROUTINE load_tensor4_from_file(s,t)
1322  CHARACTER (LEN=*), INTENT(IN) :: s
1323  TYPE(coolist4), DIMENSION(ndim), INTENT(OUT) :: t
1324  INTEGER :: i,ir,j,k,l,n,AllocStat
1325  REAL(KIND=8) :: v
1326  OPEN(30,file=s,status='old')
1327  DO i=1,ndim
1328  READ(30,*) ir,n
1329  IF (n /= 0) THEN
1330  ALLOCATE(t(i)%elems(n), stat=allocstat)
1331  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
1332  t(i)%nelems=n
1333  ENDIF
1334  DO n=1,t(i)%nelems
1335  READ(30,*) ir,j,k,l,v
1336  t(i)%elems(n)%j=j
1337  t(i)%elems(n)%k=k
1338  t(i)%elems(n)%l=l
1339  t(i)%elems(n)%v=v
1340  ENDDO
1341  END DO
1342  CLOSE(30)
1343  END SUBROUTINE load_tensor4_from_file
1344 
1345  !> Load a rank-4 tensor coolist from a file definition
1346  !> @param s Destination filename
1347  !> @param t The coolist to write
1348  SUBROUTINE write_tensor4_to_file(s,t)
1349  CHARACTER (LEN=*), INTENT(IN) :: s
1350  TYPE(coolist4), DIMENSION(ndim), INTENT(IN) :: t
1351  INTEGER :: i,j,k,l,n
1352  OPEN(30,file=s)
1353  DO i=1,ndim
1354  WRITE(30,*) i,t(i)%nelems
1355  DO n=1,t(i)%nelems
1356  j=t(i)%elems(n)%j
1357  k=t(i)%elems(n)%k
1358  l=t(i)%elems(n)%l
1359  WRITE(30,*) i,j,k,l,t(i)%elems(n)%v
1360  END DO
1361  END DO
1362  CLOSE(30)
1363  END SUBROUTINE write_tensor4_to_file
1364 
1365 END MODULE tensor
1366 
subroutine, public add_check(t, i, j, k, v, dst)
Subroutine to add element to a coolist and check for overflow. Once the t buffer tensor is full...
Definition: tensor.f90:332
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
subroutine, public add_matc_to_tensor4(i, j, src, dst)
Routine to add a matrix to a rank-4 tensor.
Definition: tensor.f90:537
Utility module.
Definition: util.f90:12
subroutine, public add_vec_ijk_to_tensor4(i, j, k, src, dst)
Routine to add a vector to a rank-4 tensor.
Definition: tensor.f90:785
subroutine, public copy_coo(src, dst)
Routine to copy a coolist.
Definition: tensor.f90:72
Statistics accumulators.
Definition: stat.f90:14
subroutine, public add_vec_jk_to_tensor(j, k, src, dst)
Routine to add a vector to a rank-3 tensor.
Definition: tensor.f90:602
subroutine, public load_tensor4_from_file(s, t)
Load a rank-4 tensor coolist from a file definition.
Definition: tensor.f90:1322
subroutine, public tensor4_to_coo4(src, dst)
Routine to convert a rank-4 tensor from matrix to coolist representation.
Definition: tensor.f90:883
subroutine, public simplify(tensor)
Routine to simplify a coolist (sparse tensor). For each index , it upper triangularize the matrix ...
Definition: tensor.f90:238
subroutine, public sparse_mul4_with_mat_jl(coolist_ijkl, mat_jl, res)
Sparse multiplication of a rank-4 tensor coolist with a matrix : .
Definition: tensor.f90:1169
subroutine, public mat_to_coo(src, dst)
Routine to convert a matrix to a tensor.
Definition: tensor.f90:94
subroutine, public scal_mul_coo(s, t)
Routine to multiply a rank-3 tensor by a scalar.
Definition: tensor.f90:1274
subroutine, public coo_to_mat_i(i, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
Definition: tensor.f90:1112
subroutine, public sparse_mul3(coolist_ijk, arr_j, arr_k, res)
Sparse multiplication of a tensor with two vectors: .
Definition: tensor.f90:129
subroutine, public coo_to_mat_j(j, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
Definition: tensor.f90:1148
Tensor utility module.
Definition: tensor.f90:18
subroutine, public tensor_to_coo(src, dst)
Routine to convert a rank-3 tensor from matrix to coolist representation.
Definition: tensor.f90:847
subroutine, public coo_to_mat_ij(src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix with i and j indices.
Definition: tensor.f90:1079
character(len=20) function, public str(k)
Convert an integer to string.
Definition: util.f90:31
subroutine, public add_vec_ikl_to_tensor4_perm(i, k, l, src, dst)
Routine to add a vector to a rank-4 tensor plus permutation.
Definition: tensor.f90:657
subroutine, public jsparse_mul(coolist_ijk, arr_j, jcoo_ij)
Sparse multiplication of two tensors to determine the Jacobian: It&#39;s implemented slightly differentl...
Definition: tensor.f90:153
subroutine, public write_tensor_to_file(s, t)
Load a rank-4 tensor coolist from a file definition.
Definition: tensor.f90:425
subroutine, public add_matc_to_tensor(i, src, dst)
Routine to add a matrix to a rank-3 tensor.
Definition: tensor.f90:474
Coordinate list. Type used to represent the sparse tensor.
Definition: tensor.f90:38
logical function, public tensor_empty(t)
Test if a rank-3 tensor coolist is empty.
Definition: tensor.f90:1288
subroutine, public sparse_mul4_with_mat_kl(coolist_ijkl, mat_kl, res)
Sparse multiplication of a rank-4 tensor coolist with a matrix : .
Definition: tensor.f90:1194
subroutine, public write_tensor4_to_file(s, t)
Load a rank-4 tensor coolist from a file definition.
Definition: tensor.f90:1349
subroutine, public jsparse_mul_mat(coolist_ijk, arr_j, jcoo_ij)
Sparse multiplication of two tensors to determine the Jacobian: It&#39;s implemented slightly differentl...
Definition: tensor.f90:196
subroutine, public load_tensor_from_file(s, t)
Load a rank-4 tensor coolist from a file definition.
Definition: tensor.f90:445
subroutine, public sparse_mul3_mat(coolist_ijk, arr_k, res)
Sparse multiplication of a rank-3 tensor coolist with a vector: . Its output is a matrix...
Definition: tensor.f90:948
subroutine, public coo_to_mat_ik(src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix with i and k indices.
Definition: tensor.f90:1063
subroutine, public add_to_tensor(src, dst)
Routine to add a rank-3 tensor to another one.
Definition: tensor.f90:350
subroutine, public print_tensor4(t)
Routine to print a rank-4 tensor coolist.
Definition: tensor.f90:922
subroutine, public sparse_mul2_j(coolist_ijk, arr_j, res)
Sparse multiplication of a 3d sparse tensor with a vectors: .
Definition: tensor.f90:1024
4d coordinate list. Type used to represent the rank-4 sparse tensor.
Definition: tensor.f90:44
subroutine, public matc_to_coo(src, dst)
Routine to convert a matrix to a rank-3 tensor.
Definition: tensor.f90:1244
subroutine, public sparse_mul3_with_mat(coolist_ijk, mat_jk, res)
Sparse multiplication of a rank-3 tensor coolist with a matrix: .
Definition: tensor.f90:1220
subroutine, public add_vec_ikl_to_tensor4(i, k, l, src, dst)
Routine to add a vector to a rank-4 tensor.
Definition: tensor.f90:726
subroutine, public sparse_mul2_k(coolist_ijk, arr_k, res)
Sparse multiplication of a rank-3 sparse tensor coolist with a vector: .
Definition: tensor.f90:1045
The model parameters module.
Definition: params.f90:18
Coordinate list element type. Elementary elements of the sparse tensors.
Definition: tensor.f90:25
subroutine, public sparse_mul2(coolist_ij, arr_j, res)
Sparse multiplication of a 2d sparse tensor with a vector: .
Definition: tensor.f90:221
logical function, public tensor4_empty(t)
Test if a rank-4 tensor coolist is empty.
Definition: tensor.f90:1304
4d coordinate list element type. Elementary elements of the 4d sparse tensors.
Definition: tensor.f90:32
subroutine, public sparse_mul4(coolist_ijkl, arr_j, arr_k, arr_l, res)
Sparse multiplication of a rank-4 tensor coolist with three vectors: .
Definition: tensor.f90:974
real(kind=8), parameter real_eps
Parameter to test the equality with zero.
Definition: tensor.f90:50
subroutine, public add_elem(t, i, j, k, v)
Subroutine to add element to a coolist.
Definition: tensor.f90:310
subroutine, public sparse_mul4_mat(coolist_ijkl, arr_k, arr_l, res)
Sparse multiplication of a tensor with two vectors: .
Definition: tensor.f90:998
subroutine, public coo_to_vec_jk(j, k, src, dst)
Routine to convert a rank-3 tensor coolist component into a vector.
Definition: tensor.f90:1129
subroutine, public print_tensor(t, s)
Routine to print a rank 3 tensor coolist.
Definition: tensor.f90:399