ROck Physics Toolbox  1.0 A microgeodynamics-based toolkit for rock physics.
global_old.f90
Go to the documentation of this file.
1 MODULE global
2  INTEGER, PARAMETER :: i4b = selected_int_kind(9)
3  INTEGER, PARAMETER :: i2b = selected_int_kind(4)
4  INTEGER, PARAMETER :: i1b = selected_int_kind(2)
5  INTEGER, PARAMETER :: dp = kind(1.0)
6  INTEGER, PARAMETER :: sp = kind(1.0d0)
7  REAL(SP), PARAMETER :: pi=3.141592653589793238462643383279502884197_sp
8  TYPE sparse
9  REAL(sp),DIMENSION(:),POINTER::a
10  INTEGER,DIMENSION(:),POINTER::ia,ja
11  INTEGER::m,n
12  END TYPE sparse
13  INTERFACE linspace
14  MODULE PROCEDURE linspace_r, linspace_i
15  END INTERFACE linspace
16 CONTAINS
17  PURE FUNCTION linspace_r(a,b,n)
18  IMPLICIT NONE
19  REAL(sp),INTENT(in)::a,b
20  INTEGER(I4B),INTENT(in)::n
21  INTEGER(I4B)::ii
22  REAL(sp),DIMENSION(n)::linspace_r
23  DO ii=1,n
24  linspace_r(ii)=a+(REAL(ii,kind=sp)-1.0_sp)/(REAL(n,kind=sp)-1.0_sp)*(b-a)
25  END DO
26  END FUNCTION linspace_r
27  PURE FUNCTION linspace_i(a,b,n)
28  IMPLICIT NONE
29  INTEGER(I4B),INTENT(in)::a,b
30  INTEGER(I4B),INTENT(in)::n
31  INTEGER(I4B)::ii
32  INTEGER(I4B),DIMENSION(n)::linspace_i
33  DO ii=1,n
34  linspace_i(ii)=a+(ii-1)*(b-a)/(n-1)
35  END DO
36  END FUNCTION linspace_i
37
38  pure function eye(n)
39  integer(I4B),intent(in)::n
40  real(sp),dimension(n,n)::eye
41  integer(I4B)::ii,jj
42  eye=0.0_sp
43  do ii=1,n
44  do jj=1,n
45  if(ii==jj)then
46  eye(ii,jj)=1.0_sp
47  end if
48  end do
49  end do
50  end function eye
51
52  FUNCTION lusol(KPACK,FPACK)
53
54  IMPLICIT NONE
55  REAL(sp),DIMENSION(:,:),INTENT(in)::kpack
56  REAL(sp),DIMENSION(:),INTENT(in)::fpack
57  REAL(sp),DIMENSION(SIZE(FPACK,1))::lusol
58  real(SP),dimension(size(KPACK,1),size(KPACK,2))::a
59  !This is a wrapper for solving a general system of linear
60  !equations using the intel mkl library
61  !KPACK*lusol=FPACK
62  INTEGER(I4B)::n,lda,info,ldb,nrhs
63  INTEGER(I4B),DIMENSION(SIZE(KPACK,1))::ipiv
64  CHARACTER(len=1)::trans
65  n=SIZE(kpack,1)
66  lda=n
67  lusol=fpack
68  a=kpack
69  !LU factorize the matrix first
70  CALL dgetrf(n,n,a,lda,ipiv,info)
71  trans='N'
72  nrhs=1
73  ldb=n
74  !Then solve
75  CALL dgetrs(trans,n,nrhs,a,lda,ipiv,lusol,ldb,info)
76  print*,' lusol:error',sum(abs(matmul(kpack,lusol)-fpack))
77
78  END FUNCTION lusol
79
80  FUNCTION tridag(a,b,c,r)
81  IMPLICIT NONE
82  REAL(sp),DIMENSION(:),INTENT(in)::a,b,c,r
83  REAL(sp),DIMENSION(SIZE(b,1)-1)::du2
84  REAL(sp),DIMENSION(SIZE(r,1))::tridag,sol
85  INTEGER(I4B)::n,lda,info,ldb,nrhs
86  INTEGER(I4B),DIMENSION(SIZE(b,1))::ipiv
87  CHARACTER(len=1)::trans
88  n=SIZE(b,1)
89  sol=r
90  !LU factorize the matrix first
91  CALL dgttrf(n,a,b,c,du2,ipiv,info)
92  trans='N'
93  nrhs=1
94  ldb=n
95  !Then solve
96
97  CALL dgttrs(trans,n,nrhs,a,b,c,du2,ipiv,sol,ldb,info)
98  tridag=sol
99  END FUNCTION tridag
100
101  FUNCTION gmres(KPACK,FPACK)
102
103  IMPLICIT NONE
104  REAL(sp),DIMENSION(:,:),INTENT(in)::kpack
105  REAL(sp),DIMENSION(:),INTENT(in)::fpack
106  REAL(sp),DIMENSION(SIZE(FPACK,1))::gmres,rhs,soln
107  type(sparse)::sprs
108  integer::siz
109  siz=128
110  rhs=fpack
111  soln=0.0
112  sprs=dnscsr(kpack,1.0e-6_sp)
113  call fgmres_nopc(sprs%IA,sprs%JA,sprs%A,rhs,soln,SIZE(fpack,1),siz)
114  gmres=soln
115  end FUNCTION gmres
116
117  SUBROUTINE fgmres_nopc(IA,JA,A,RHS,COMPUTED_SOLUTION,N,SIZ)
118
119  INTEGER,DIMENSION(:),INTENT(inout)::ia,ja
120  REAL(sp),DIMENSION(:),INTENT(inout)::a,computed_solution,rhs
121  INTEGER,INTENT(in)::n,siz
122  INTEGER, PARAMETER::expected_itercount=5
123
124  INTEGER::ipar(siz),itercount,rci_request, i
125  REAL(sp)::dpar(siz),tmp(n*(2*n+1)+(n*(n+9))/2+1)
126  REAL(sp):: dvar,expected_solution(n)
127
128  !---------------------------------------------------------------------------
129  ! Initialize variables and the right hand side through matrix-vector product
130  !---------------------------------------------------------------------------
131  CALL mkl_dcsrgemv('N', n, a, ia, ja, expected_solution, rhs)
132
133  !---------------------------------------------------------------------------
134  ! Initialize the initial guess
135  !---------------------------------------------------------------------------
136  DO i=1,n
137  computed_solution(i)=1.0_sp
138  ENDDO
139  !---------------------------------------------------------------------------
140  ! Initialize the solver
141  !---------------------------------------------------------------------------
142  CALL dfgmres_init(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
143  !IF (RCI_REQUEST.NE.0) GOTO 999
144  !---------------------------------------------------------------------------
145  ! Set the desired parameters:
146  ! LOGICAL parameters:
147  ! do residual stopping test
148  ! do not request for the user defined stopping test
149  ! do the check of the norm of the next generated vector automatically
150  ! DOUBLE PRECISION parameters
151  ! set the relative tolerance to 1.0D-3 instead of default value 1.0D-6
152  !---------------------------------------------------------------------------
153  ipar(9)=1
154  ipar(10)=0
155  ipar(12)=1
156  dpar(1)=1.0e-3_sp
157  !---------------------------------------------------------------------------
158  ! Check the correctness and consistency of the newly set parameters
159  !---------------------------------------------------------------------------
160  CALL dfgmres_check(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
161
162  ! Reverse Communication starts here
163  !---------------------------------------------------------------------------
164 1 CALL dfgmres(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
165  !---------------------------------------------------------------------------
166  ! If RCI_REQUEST=0, then the solution was found with the required precision
167  !---------------------------------------------------------------------------
168  IF (rci_request.EQ.0) goto 3
169  !---------------------------------------------------------------------------
170  ! If RCI_REQUEST=1, then compute the vector A*TMP(IPAR(22))
171  ! and put the result in vector TMP(IPAR(23))
172  !---------------------------------------------------------------------------
173  IF (rci_request.EQ.1) THEN
174  CALL mkl_dcsrgemv('N',n, a, ia, ja, tmp(ipar(22)), tmp(ipar(23)))
175  goto 1
176  !---------------------------------------------------------------------------
177  ! If RCI_REQUEST=anything else, then DFGMRES subroutine failed
178  ! to compute the solution vector: COMPUTED_SOLUTION(N)
179  !---------------------------------------------------------------------------
180  ELSE
181  goto 999
182  ENDIF
183  !---------------------------------------------------------------------------
184  ! Reverse Communication ends here
185  ! Get the current iteration number and the FGMRES solution (DO NOT FORGET to
186  ! call DFGMRES_GET routine as COMPUTED_SOLUTION is still containing
187  ! the initial guess!)
188  !---------------------------------------------------------------------------
189 3 CALL dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)
190  !---------------------------------------------------------------------------
191  !---------------------------------------------------------------------------
192  ! Print solution vector: COMPUTED_SOLUTION(N) and
193  ! the number of iterations: ITERCOUNT
194  !---------------------------------------------------------------------------
195
196
197  !---------------------------------------------------------------------------
198  ! Release internal MKL memory that might be used for computations
199  ! NOTE: It is important to call the routine below to avoid memory leaks
200  ! unless you disable MKL Memory Manager
201  !---------------------------------------------------------------------------
202  CALL mkl_freebuffers
203
204  dvar=dnrm2(n,expected_solution,1)
205
206  !---------------------------------------------------------------------------
207  ! Release internal MKL memory that might be used for computations
208  ! NOTE: It is important to call the routine below to avoid memory leaks
209  ! unless you disable MKL Memory Manager
210  !---------------------------------------------------------------------------
211 999 WRITE( *,'(A,A,I5)') 'This example FAILED as the solver has',&
212  & ' returned the ERROR code', rci_request
213  CALL mkl_freebuffers
214  ! STOP 1
215  END SUBROUTINE fgmres_nopc
216
217 FUNCTION dnscsr(inmat,tolin)
218 REAL(sp),DIMENSION(:,:),INTENT(in)::inmat
219 real(sp),optional::tolin
220 TYPE(sparse)::dnscsr
221 integer::nz,m,n,job(6),info
222 REAL(sp)::tol
223 REAL(sp),dimension(size(inmat,1),size(inmat,2))::mat
224 !! set the tolerance for zero search
225 IF(present(tolin)) THEN
226  tol=tolin
227  ELSE
228  tol=1.0e-8_sp
229  END IF
230 mat=inmat
231 !Zero out the entries that are below tolerance level
232 where(abs(inmat)<tol)
233 mat=0.0_sp
234 end where
235 ! Count the number of nonzero entries in the incoming matrix
236 nz=sparscalc(mat,tol)
237 !Number of rows in the incoming matrix
238 m=size(mat,1)
239 !Number of columns in the incoming matrix
240 n=size(mat,2)
241 ! Allocate the arrays for sparse storage
242 ALLOCATE(dnscsr%A(nz),dnscsr%IA(m+1),dnscsr%JA(nz))
243
244 ! Set the job preference
245 job(1)=0 ! Convert regular matrix to CSR
246 job(2)=1 ! Use one-based matrix notation
247 job(3)=1 ! Use one-based indexing in CSR format
248 job(4)=2 ! Input the whole matrix
249 job(5)=nz! Maximum number of nonzero elements
250 job(6)=1 ! Return Ia, Ja, and AA
251
252 call mkl_ddnscsr(job,m,n,mat,m,dnscsr%A,dnscsr%JA,dnscsr%IA,info)
253
254 END FUNCTION dnscsr
255
256 FUNCTION sparscalc(mat,tolin)
257  ! This function calculates the number of nonzero entries in
258  ! a sparse matrix. By default, numbers smaller than 1.0e-8
259  ! are considered zero, but the default can be overriden by
260  ! supplying the optional second paramter
261
262  implicit none
263  REAL(sp),DIMENSION(:,:),INTENT(in)::mat
264  REAL(sp),OPTIONAL::tolin
265  REAL(sp)::tol
266  INTEGER::nn,sparscalc
267  INTEGER,DIMENSION(SIZE(mat,1),SIZE(mat,2))::tmp
268  tmp=0
269  IF(present(tolin)) THEN
270  tol=tolin
271  ELSE
272  tol=1.0e-8_sp
273  END IF
274
275  WHERE(abs(mat)>tol)
276  tmp=1
277  END WHERE
278  sparscalc=sum(tmp)
279 END FUNCTION sparscalc
280
281  FUNCTION filename(f,str)
282  INTEGER(I4B),INTENT(in)::f
283  character(len=3),intent(in)::str
284  CHARACTER(len=40)::frmt
285  CHARACTER(len=18)::filename
286
287  IF(f<10) THEN
288  frmt='("../data/",a3,"00",i1,".vtk")'
289  WRITE(filename,frmt),str,f
290  elseif((f<100).AND.(f>=10)) THEN
291  frmt='("../data/",a3,"0",i2,".vtk")'
292  WRITE(filename,frmt),str,f
293  elseif((f<1000).AND.(f>=100)) THEN
294  frmt='("../data/",a3,i3,".vtk")'
295  WRITE(filename,frmt),str,f
296  END IF
297  END FUNCTION filename
298
299
300 END MODULE global