FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 
25 contains
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 
149 end module hecmw_precond_diag_nn
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
subroutine, public hecmw_precond_diag_nn_setup(hecMAT)
subroutine, public hecmw_precond_diag_nn_clear()
subroutine, public hecmw_precond_diag_nn_apply(WW, NDOF)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal