50 REAL(KIND=8),
PARAMETER ::
real_eps = 2.2204460492503131e-16
72 TYPE(
coolist),
DIMENSION(ndim),
INTENT(IN) :: src
73 TYPE(
coolist),
DIMENSION(ndim),
INTENT(OUT) :: dst
74 INTEGER :: i,j,AllocStat
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 ! ***" 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
85 dst(i)%nelems=src(i)%nelems
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
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 ! ***" 112 dst(i)%elems(n)%v=src(i,j)
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
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)
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
157 INTEGER :: i,j,k,n,nj,AllocStat
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 ! ***" 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
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)
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)
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
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)
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
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)
238 TYPE(
coolist),
DIMENSION(ndim),
INTENT(INOUT):: tensor
240 INTEGER :: li,lii,liii,n
248 &%elems(lii)%k)).OR.((j==
tensor(i)%elems(lii)%k).AND.(k==
tensor(i)%elems(lii)%j)))
THEN 281 if (li >
tensor(i)%nelems)
THEN 310 TYPE(
coolist),
DIMENSION(ndim),
INTENT(INOUT) :: t
311 INTEGER,
INTENT(IN) :: i,j,k
312 REAL(KIND=8),
INTENT(IN) :: v
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
338 IF (t(i)%nelems==
size(t(i)%elems))
THEN 350 TYPE(
coolist),
DIMENSION(ndim),
INTENT(IN) :: src
351 TYPE(
coolist),
DIMENSION(ndim),
INTENT(INOUT) :: dst
353 INTEGER :: i,j,n,AllocStat
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 ! ***" 362 ALLOCATE(dst(i)%elems(src(i)%nelems),
stat=allocstat)
363 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 367 ALLOCATE(celems(n),
stat=allocstat)
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
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 ! ***" 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
382 DEALLOCATE(celems,
stat=allocstat)
383 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 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
390 dst(i)%nelems=src(i)%nelems+n
400 TYPE(
coolist),
DIMENSION(ndim),
INTENT(IN) :: t
401 CHARACTER,
INTENT(IN),
OPTIONAL :: s
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
425 CHARACTER (LEN=*),
INTENT(IN) :: s
426 TYPE(
coolist),
DIMENSION(ndim),
INTENT(IN) :: t
430 WRITE(30,*) i,t(i)%nelems
434 WRITE(30,*) i,j,k,t(i)%elems(n)%v
445 CHARACTER (LEN=*),
INTENT(IN) :: s
446 TYPE(
coolist),
DIMENSION(ndim),
INTENT(OUT) :: t
447 INTEGER :: i,ir,j,k,n,AllocStat
449 OPEN(30,file=s,status=
'old')
453 ALLOCATE(t(i)%elems(n),
stat=allocstat)
454 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 474 INTEGER,
INTENT(IN) :: i
475 REAL(KIND=8),
DIMENSION(ndim,ndim),
INTENT(IN) :: src
476 TYPE(
coolist),
DIMENSION(ndim),
INTENT(INOUT) :: dst
478 INTEGER :: j,k,r,n,nsrc,AllocStat
483 IF (abs(src(j,k))>
real_eps) nsrc=nsrc+1
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 ! ***" 492 ALLOCATE(dst(i)%elems(nsrc),
stat=allocstat)
493 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 497 ALLOCATE(celems(n),
stat=allocstat)
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
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 ! ***" 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
512 DEALLOCATE(celems,
stat=allocstat)
513 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 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)
537 INTEGER,
INTENT(IN) :: i,j
538 REAL(KIND=8),
DIMENSION(ndim,ndim),
INTENT(IN) :: src
539 TYPE(
coolist4),
DIMENSION(ndim),
INTENT(INOUT) :: dst
541 INTEGER :: k,l,r,n,nsrc,AllocStat
546 IF (abs(src(k,l))>
real_eps) nsrc=nsrc+1
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 ! ***" 555 ALLOCATE(dst(i)%elems(nsrc),
stat=allocstat)
556 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 560 ALLOCATE(celems(n),
stat=allocstat)
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
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 ! ***" 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
577 DEALLOCATE(celems,
stat=allocstat)
578 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 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)
602 INTEGER,
INTENT(IN) :: j,k
603 REAL(KIND=8),
DIMENSION(ndim),
INTENT(IN) :: src
604 TYPE(
coolist),
DIMENSION(ndim),
INTENT(INOUT) :: dst
606 INTEGER :: i,l,r,n,nsrc,AllocStat
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 ! ***" 616 ALLOCATE(dst(i)%elems(nsrc),
stat=allocstat)
617 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 621 ALLOCATE(celems(n),
stat=allocstat)
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
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 ! ***" 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
636 DEALLOCATE(celems,
stat=allocstat)
637 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 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)
657 INTEGER,
INTENT(IN) :: i,k,l
658 REAL(KIND=8),
DIMENSION(ndim),
INTENT(IN) :: src
659 TYPE(
coolist4),
DIMENSION(ndim),
INTENT(INOUT) :: dst
661 INTEGER :: j,ne,r,n,nsrc,AllocStat
665 IF (abs(src(j))>
real_eps) nsrc=nsrc+1
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 ! ***" 673 ALLOCATE(dst(i)%elems(nsrc),
stat=allocstat)
674 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 678 ALLOCATE(celems(n),
stat=allocstat)
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
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 ! ***" 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
695 DEALLOCATE(celems,
stat=allocstat)
696 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 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)
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)
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)
726 INTEGER,
INTENT(IN) :: i,k,l
727 REAL(KIND=8),
DIMENSION(ndim),
INTENT(IN) :: src
728 TYPE(
coolist4),
DIMENSION(ndim),
INTENT(INOUT) :: dst
730 INTEGER :: j,ne,r,n,nsrc,AllocStat
734 IF (abs(src(j))>
real_eps) nsrc=nsrc+1
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 ! ***" 742 ALLOCATE(dst(i)%elems(nsrc),
stat=allocstat)
743 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 747 ALLOCATE(celems(n),
stat=allocstat)
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
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 ! ***" 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
764 DEALLOCATE(celems,
stat=allocstat)
765 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 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)
785 INTEGER,
INTENT(IN) :: i,j,k
786 REAL(KIND=8),
DIMENSION(ndim),
INTENT(IN) :: src
787 TYPE(
coolist4),
DIMENSION(ndim),
INTENT(INOUT) :: dst
789 INTEGER :: l,ne,r,n,nsrc,AllocStat
793 IF (abs(src(l))>
real_eps) nsrc=nsrc+1
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 ! ***" 801 ALLOCATE(dst(i)%elems(nsrc),
stat=allocstat)
802 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 806 ALLOCATE(celems(n),
stat=allocstat)
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
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 ! ***" 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
823 DEALLOCATE(celems,
stat=allocstat)
824 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 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)
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
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 ! ***" 869 dst(i)%elems(n)%v=src(i,j,k)
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
892 IF (abs(src(i,j,k,l))>
real_eps) n=n+1
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 ! ***" 904 IF (abs(src(i,j,k,l))>
real_eps)
THEN 909 dst(i)%elems(n)%v=src(i,j,k,l)
923 TYPE(
coolist4),
DIMENSION(ndim),
INTENT(IN) :: t
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
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
954 DO n=1,coolist_ijk(i)%nelems
955 j=coolist_ijk(i)%elems(n)%j
957 k=coolist_ijk(i)%elems(n)%k
958 res(i,j) = res(i,j) + coolist_ijk(i)%elems(n)%v * arr_k(k)
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
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)
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
1004 DO n=1,coolist_ijkl(i)%nelems
1005 j=coolist_ijkl(i)%elems(n)%j
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)
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
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)
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
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)
1064 REAL(KIND=8),
DIMENSION(ndim,ndim),
INTENT(OUT) :: dst
1069 DO n=1,src(i)%nelems
1070 dst(i,src(i)%elems(n)%k)=src(i)%elems(n)%v
1080 REAL(KIND=8),
DIMENSION(ndim,ndim),
INTENT(OUT) :: dst
1085 DO n=1,src(i)%nelems
1086 dst(i,src(i)%elems(n)%j)=src(i)%elems(n)%v
1112 INTEGER,
INTENT(IN) :: i
1113 TYPE(
coolist),
DIMENSION(ndim),
INTENT(IN) :: src
1114 REAL(KIND=8),
DIMENSION(ndim,ndim),
INTENT(OUT) :: dst
1118 DO n=1,src(i)%nelems
1119 dst(src(i)%elems(n)%j,src(i)%elems(n)%k)=src(i)%elems(n)%v
1129 INTEGER,
INTENT(IN) :: j,k
1130 TYPE(
coolist),
DIMENSION(ndim),
INTENT(IN) :: src
1131 REAL(KIND=8),
DIMENSION(ndim),
INTENT(OUT) :: dst
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
1148 INTEGER,
INTENT(IN) :: j
1149 TYPE(
coolist),
DIMENSION(ndim),
INTENT(IN) :: src
1150 REAL(KIND=8),
DIMENSION(ndim,ndim),
INTENT(OUT) :: dst
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
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
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
1181 res(i,k) = res(i,k) + coolist_ijkl(i)%elems(n)%v * mat_jl(j,l)
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
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
1206 res(i,j) = res(i,j) + coolist_ijkl(i)%elems(n)%v * mat_kl(k,l)
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
1227 DO n=1,coolist_ijk(i)%nelems
1228 j=coolist_ijk(i)%elems(n)%j
1229 k=coolist_ijk(i)%elems(n)%k
1231 res(i) = res(i) + coolist_ijk(i)%elems(n)%v * mat_jk(j,k)
1244 REAL(KIND=8),
DIMENSION(ndim,ndim),
INTENT(IN) :: src
1245 TYPE(
coolist),
DIMENSION(ndim),
INTENT(OUT) :: dst
1246 INTEGER :: i,j,n,AllocStat
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 ! ***" 1262 dst(i)%elems(n)%v=src(i,j)
1274 REAL(KIND=8),
INTENT(IN) :: s
1275 TYPE(
coolist),
DIMENSION(ndim),
INTENT(INOUT) :: t
1280 t(i)%elems(li)%v=s*t(i)%elems(li)%v
1289 LOGICAL :: tensor_empty
1293 IF (t(i)%nelems /= 0)
THEN 1294 tensor_empty=.false.
1305 LOGICAL :: tensor4_empty
1307 tensor4_empty=.true.
1309 IF (t(i)%nelems /= 0)
THEN 1310 tensor4_empty=.false.
1322 CHARACTER (LEN=*),
INTENT(IN) :: s
1323 TYPE(
coolist4),
DIMENSION(ndim),
INTENT(OUT) :: t
1324 INTEGER :: i,ir,j,k,l,n,AllocStat
1326 OPEN(30,file=s,status=
'old')
1330 ALLOCATE(t(i)%elems(n),
stat=allocstat)
1331 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 1335 READ(30,*) ir,j,k,l,v
1349 CHARACTER (LEN=*),
INTENT(IN) :: s
1350 TYPE(
coolist4),
DIMENSION(ndim),
INTENT(IN) :: t
1351 INTEGER :: i,j,k,l,n
1354 WRITE(30,*) i,t(i)%nelems
1359 WRITE(30,*) i,j,k,l,t(i)%elems(n)%v
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...
integer ndim
Number of variables (dimension of the model)
subroutine, public add_matc_to_tensor4(i, j, src, dst)
Routine to add a matrix to a rank-4 tensor.
subroutine, public add_vec_ijk_to_tensor4(i, j, k, src, dst)
Routine to add a vector to a rank-4 tensor.
subroutine, public copy_coo(src, dst)
Routine to copy a coolist.
subroutine, public add_vec_jk_to_tensor(j, k, src, dst)
Routine to add a vector to a rank-3 tensor.
subroutine, public load_tensor4_from_file(s, t)
Load a rank-4 tensor coolist from a file definition.
subroutine, public tensor4_to_coo4(src, dst)
Routine to convert a rank-4 tensor from matrix to coolist representation.
subroutine, public simplify(tensor)
Routine to simplify a coolist (sparse tensor). For each index , it upper triangularize the matrix ...
subroutine, public sparse_mul4_with_mat_jl(coolist_ijkl, mat_jl, res)
Sparse multiplication of a rank-4 tensor coolist with a matrix : .
subroutine, public mat_to_coo(src, dst)
Routine to convert a matrix to a tensor.
subroutine, public scal_mul_coo(s, t)
Routine to multiply a rank-3 tensor by a scalar.
subroutine, public coo_to_mat_i(i, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
subroutine, public sparse_mul3(coolist_ijk, arr_j, arr_k, res)
Sparse multiplication of a tensor with two vectors: .
subroutine, public coo_to_mat_j(j, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
subroutine, public tensor_to_coo(src, dst)
Routine to convert a rank-3 tensor from matrix to coolist representation.
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.
character(len=20) function, public str(k)
Convert an integer to string.
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.
subroutine, public jsparse_mul(coolist_ijk, arr_j, jcoo_ij)
Sparse multiplication of two tensors to determine the Jacobian: It's implemented slightly differentl...
subroutine, public write_tensor_to_file(s, t)
Load a rank-4 tensor coolist from a file definition.
subroutine, public add_matc_to_tensor(i, src, dst)
Routine to add a matrix to a rank-3 tensor.
Coordinate list. Type used to represent the sparse tensor.
logical function, public tensor_empty(t)
Test if a rank-3 tensor coolist is empty.
subroutine, public sparse_mul4_with_mat_kl(coolist_ijkl, mat_kl, res)
Sparse multiplication of a rank-4 tensor coolist with a matrix : .
subroutine, public write_tensor4_to_file(s, t)
Load a rank-4 tensor coolist from a file definition.
subroutine, public jsparse_mul_mat(coolist_ijk, arr_j, jcoo_ij)
Sparse multiplication of two tensors to determine the Jacobian: It's implemented slightly differentl...
subroutine, public load_tensor_from_file(s, t)
Load a rank-4 tensor coolist from a file definition.
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...
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.
subroutine, public add_to_tensor(src, dst)
Routine to add a rank-3 tensor to another one.
subroutine, public print_tensor4(t)
Routine to print a rank-4 tensor coolist.
subroutine, public sparse_mul2_j(coolist_ijk, arr_j, res)
Sparse multiplication of a 3d sparse tensor with a vectors: .
4d coordinate list. Type used to represent the rank-4 sparse tensor.
subroutine, public matc_to_coo(src, dst)
Routine to convert a matrix to a rank-3 tensor.
subroutine, public sparse_mul3_with_mat(coolist_ijk, mat_jk, res)
Sparse multiplication of a rank-3 tensor coolist with a matrix: .
subroutine, public add_vec_ikl_to_tensor4(i, k, l, src, dst)
Routine to add a vector to a rank-4 tensor.
subroutine, public sparse_mul2_k(coolist_ijk, arr_k, res)
Sparse multiplication of a rank-3 sparse tensor coolist with a vector: .
The model parameters module.
Coordinate list element type. Elementary elements of the sparse tensors.
subroutine, public sparse_mul2(coolist_ij, arr_j, res)
Sparse multiplication of a 2d sparse tensor with a vector: .
logical function, public tensor4_empty(t)
Test if a rank-4 tensor coolist is empty.
4d coordinate list element type. Elementary elements of the 4d sparse tensors.
subroutine, public sparse_mul4(coolist_ijkl, arr_j, arr_k, arr_l, res)
Sparse multiplication of a rank-4 tensor coolist with three vectors: .
real(kind=8), parameter real_eps
Parameter to test the equality with zero.
subroutine, public add_elem(t, i, j, k, v)
Subroutine to add element to a coolist.
subroutine, public sparse_mul4_mat(coolist_ijkl, arr_k, arr_l, res)
Sparse multiplication of a tensor with two vectors: .
subroutine, public coo_to_vec_jk(j, k, src, dst)
Routine to convert a rank-3 tensor coolist component into a vector.
subroutine, public print_tensor(t, s)
Routine to print a rank 3 tensor coolist.