FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_precond_DIAG_nn.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
6!C
7!C***
8!C*** module hecmw_precond_DIAG_nn
9!C***
10!C
12 use hecmw_util
13
14 private
15
19
20 integer(kind=kint) :: N
21 real(kind=kreal), pointer :: alu(:) => null()
22
23 logical, save :: INITIALIZED = .false.
24
25contains
26
27 subroutine hecmw_precond_diag_nn_setup(hecMAT)
29 implicit none
30 type(hecmwst_matrix), intent(inout) :: hecmat
31 integer(kind=kint ) :: np,ndof,ndof2
32 real (kind=kreal) :: sigma_diag
33 real(kind=kreal), pointer:: d(:)
34
35 real (kind=kreal):: alutmp(hecmat%NDOF,hecmat%NDOF), pw(hecmat%NDOF)
36 integer(kind=kint ):: i, j, k, ii
37
38 if (initialized) then
39 if (hecmat%Iarray(98) == 0 .and. hecmat%Iarray(97) == 0) return
41 endif
42
43 n = hecmat%N
44 ndof = hecmat%NDOF
45 ndof2 = ndof*ndof
46 np = hecmat%NP
47 d => hecmat%D
48 sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
49
50 allocate(alu(ndof2*np))
51 alu = 0.d0
52
53 do ii= 1, n
54 do i = 1, ndof2
55 alu(ndof2*(ii-1)+i)=d(ndof2*(ii-1)+i)
56 end do
57 enddo
58
59 if (hecmat%cmat%n_val.gt.0) then
60 do k= 1, hecmat%cmat%n_val
61 if (hecmat%cmat%pair(k)%i.ne.hecmat%cmat%pair(k)%j) cycle
62 ii = hecmat%cmat%pair(k)%i
63 do i = 1,ndof
64 do j = 1,ndof
65 alu(ndof2*(ii-1)+(i-1)*ndof+j) = alu(ndof2*(ii-1)+(i-1)*ndof+j) + hecmat%cmat%pair(k)%val(i, j)
66 enddo
67 enddo
68 enddo
69 endif
70
71 !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,NDOF,NDOF2,ALU,SIGMA_DIAG)
72 !$omp do
73 do ii= 1, n
74
75 do i = 1, ndof
76 do j = 1, ndof
77 alutmp(i,j) = alu(ndof2*(ii-1)+(i-1)*ndof+j)
78 if (i==j) alutmp(i,j)=alutmp(i,j)*sigma_diag
79 end do
80 end do
81 do k= 1, ndof
82 alutmp(k,k)= 1.d0/alutmp(k,k)
83 do i= k+1, ndof
84 alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
85 do j= k+1, ndof
86 pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
87 enddo
88 do j= k+1, ndof
89 alutmp(i,j)= pw(j)
90 enddo
91 enddo
92 enddo
93 do i = 1, ndof
94 do j = 1, ndof
95 alu(ndof2*(ii-1)+(i-1)*ndof+j)= alutmp(i,j)
96 end do
97 end do
98 enddo
99 !$omp end do
100 !$omp end parallel
101 initialized = .true.
102 hecmat%Iarray(98) = 0 ! symbolic setup done
103 hecmat%Iarray(97) = 0 ! numerical setup done
104
105 end subroutine hecmw_precond_diag_nn_setup
106
107 subroutine hecmw_precond_diag_nn_apply(WW,NDOF)
108 implicit none
109 real(kind=kreal), intent(inout) :: ww(:)
110 integer(kind=kint) :: i,j,k,ndof
111 real(kind=kreal) :: x(ndof)
112
113 !C
114 !C== Block SCALING
115 !$omp parallel default(none),private(i,j,k,X),shared(N,WW,ALU,NDOF)
116 !$omp do
117 do i= 1, n
118 do j=1,ndof
119 x(j)=ww(ndof*(i-1)+j)
120 end do
121 do j=2,ndof
122 do k = 1,j-1
123 x(j)=x(j)-alu(ndof*ndof*(i-1)+ndof*(j-1)+k )*x(k)
124 end do
125 end do
126 do j=ndof,1,-1
127 do k = ndof,j+1,-1
128 x(j)=x(j)-alu(ndof*ndof*(i-1)+ndof*(j-1)+k )*x(k)
129 end do
130 x(j)=alu(ndof*ndof*(i-1)+(ndof+1)*(j-1)+1 )*x(j)
131 end do
132 do j=1,ndof
133 ww(ndof*(i-1)+j)=x(j)
134 end do
135
136 enddo
137 !$omp end do
138 !$omp end parallel
139
140 end subroutine hecmw_precond_diag_nn_apply
141
143 implicit none
144 if (associated(alu)) deallocate(alu)
145 nullify(alu)
146 initialized = .false.
147 end subroutine hecmw_precond_diag_nn_clear
148
149end module hecmw_precond_diag_nn
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecmat)
subroutine, public hecmw_precond_diag_nn_clear()
subroutine, public hecmw_precond_diag_nn_setup(hecmat)
subroutine, public hecmw_precond_diag_nn_apply(ww, ndof)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal