ROck Physics Toolbox  1.0 A microgeodynamics-based toolkit for rock physics.
global.f90
Go to the documentation of this file.
1 MODULE global
2  !USE LAPACK95
3  INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
4  INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
5  INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
6  INTEGER, PARAMETER :: DP = KIND(1.0)
7  INTEGER, PARAMETER :: SP = KIND(1.0D0)
8  REAL(SP), PARAMETER :: PI=3.141592653589793238462643383279502884197_sp
9  TYPE sparse
10  REAL(sp),DIMENSION(:),POINTER::A
11  INTEGER,DIMENSION(:),POINTER::IA,JA
12  INTEGER::m,n
13  END TYPE sparse
14  INTERFACE linspace
15  MODULE PROCEDURE linspace_r, linspace_i
16  END INTERFACE linspace
17
18 CONTAINS
19  PURE FUNCTION linspace_r(a,b,n)
20  IMPLICIT NONE
21  REAL(sp),INTENT(in)::a,b
22  INTEGER(I4B),INTENT(in)::n
23  INTEGER(I4B)::ii
24  REAL(sp),DIMENSION(n)::linspace_r
25  DO ii=1,n
26  linspace_r(ii)=a+(REAL(ii,kind=sp)-1.0_sp)/(REAL(n,kind=sp)-1.0_sp)*(b-a)
27  END DO
28  END FUNCTION linspace_r
29  PURE FUNCTION linspace_i(a,b,n)
30  IMPLICIT NONE
31  INTEGER(I4B),INTENT(in)::a,b
32  INTEGER(I4B),INTENT(in)::n
33  INTEGER(I4B)::ii
34  INTEGER(I4B),DIMENSION(n)::linspace_i
35  DO ii=1,n
36  linspace_i(ii)=a+(ii-1)*(b-a)/(n-1)
37  END DO
38  END FUNCTION linspace_i
39
40  PURE FUNCTION eye(n)
41  INTEGER(I4B),INTENT(in)::n
42  REAL(sp),DIMENSION(n,n)::eye
43  INTEGER(I4B)::ii,jj
44  eye=0.0_sp
45  DO ii=1,n
46  DO jj=1,n
47  IF(ii==jj)THEN
48  eye(ii,jj)=1.0_sp
49  END IF
50  END DO
51  END DO
52  END FUNCTION eye
53
54  FUNCTION lusol(KPACK,FPACK)
55
56  IMPLICIT NONE
57  REAL(sp),DIMENSION(:,:),INTENT(in)::kpack
58  REAL(sp),DIMENSION(:),INTENT(in)::fpack
59  REAL(sp),DIMENSION(SIZE(FPACK,1))::lusol
60  REAL(SP),DIMENSION(SIZE(KPACK,1),SIZE(KPACK,2))::a
61  !This is a wrapper for solving a general system of linear
62  !equations using the intel mkl library
63  !KPACK*lusol=FPACK
64  INTEGER(I4B)::n,lda,info,ldb,nrhs
65  INTEGER(I4B),DIMENSION(SIZE(KPACK,1))::ipiv
66  CHARACTER(len=1)::trans
67  n=SIZE(kpack,1)
68  lda=n
69  lusol=fpack
70  a=kpack
71  !LU factorize the matrix first
72  CALL dgetrf(n,n,a,lda,ipiv,info)
73  trans='N'
74  nrhs=1
75  ldb=n
76  !Then solve
77  CALL dgetrs(trans,n,nrhs,a,lda,ipiv,lusol,ldb,info)
78  print*,' lusol:error',sum(abs(matmul(kpack,lusol)-fpack))
79
80  END FUNCTION lusol
81
82
83  FUNCTION tridag(a,b,c,r)
84  IMPLICIT NONE
85  REAL(sp),DIMENSION(:),INTENT(in)::a,b,c,r
86  REAL(sp),DIMENSION(SIZE(b,1)-1)::du2
87  REAL(sp),DIMENSION(SIZE(r,1))::tridag,sol
88  INTEGER(I4B)::n,lda,info,ldb,nrhs
89  INTEGER(I4B),DIMENSION(SIZE(b,1))::ipiv
90  CHARACTER(len=1)::trans
91  n=SIZE(b,1)
92  sol=r
93  !LU factorize the matrix first
94  CALL dgttrf(n,a,b,c,du2,ipiv,info)
95  trans='N'
96  nrhs=1
97  ldb=n
98  !Then solve
99
100  CALL dgttrs(trans,n,nrhs,a,b,c,du2,ipiv,sol,ldb,info)
101  tridag=sol
102  END FUNCTION tridag
103
104  FUNCTION gmres(KPACK,FPACK)
105
106  IMPLICIT NONE
107  REAL(sp),DIMENSION(:,:),INTENT(in)::kpack
108  REAL(sp),DIMENSION(:),INTENT(in)::fpack
109  REAL(sp),DIMENSION(SIZE(FPACK,1))::gmres,rhs,soln
110  TYPE(sparse)::sprs
111  INTEGER::siz
112  siz=128
113  rhs=fpack
114  soln=0.0
115  sprs=dnscsr(kpack,1.0e-6_sp)
116  CALL fgmres_nopc(sprs%IA,sprs%JA,sprs%A,rhs,soln,SIZE(fpack,1),siz)
117  gmres=soln
118  END FUNCTION gmres
119
120  SUBROUTINE fgmres_nopc(IA,JA,A,RHS,COMPUTED_SOLUTION,N,SIZ)
121
122  INTEGER,DIMENSION(:),INTENT(inout)::ia,ja
123  REAL(sp),DIMENSION(:),INTENT(inout)::a,computed_solution,rhs
124  INTEGER,INTENT(in)::n,siz
125  INTEGER, PARAMETER::expected_itercount=5
126
127  INTEGER::ipar(siz),itercount,rci_request, i
128  REAL(sp)::dpar(siz),tmp(n*(2*n+1)+(n*(n+9))/2+1)
129  REAL(sp):: dvar,expected_solution(n)
130
131  !---------------------------------------------------------------------------
132  ! Initialize variables and the right hand side through matrix-vector product
133  !---------------------------------------------------------------------------
134  CALL mkl_dcsrgemv('N', n, a, ia, ja, expected_solution, rhs)
135
136  !---------------------------------------------------------------------------
137  ! Initialize the initial guess
138  !---------------------------------------------------------------------------
139  DO i=1,n
140  computed_solution(i)=1.0_sp
141  ENDDO
142  !---------------------------------------------------------------------------
143  ! Initialize the solver
144  !---------------------------------------------------------------------------
145  CALL dfgmres_init(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
146  !IF (RCI_REQUEST.NE.0) GOTO 999
147  !---------------------------------------------------------------------------
148  ! Set the desired parameters:
149  ! LOGICAL parameters:
150  ! do residual stopping test
151  ! do not request for the user defined stopping test
152  ! do the check of the norm of the next generated vector automatically
153  ! DOUBLE PRECISION parameters
154  ! set the relative tolerance to 1.0D-3 instead of default value 1.0D-6
155  !---------------------------------------------------------------------------
156  ipar(9)=1
157  ipar(10)=0
158  ipar(12)=1
159  dpar(1)=1.0e-3_sp
160  !---------------------------------------------------------------------------
161  ! Check the correctness and consistency of the newly set parameters
162  !---------------------------------------------------------------------------
163  CALL dfgmres_check(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
164
165  ! Reverse Communication starts here
166  !---------------------------------------------------------------------------
167 1 CALL dfgmres(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
168  !---------------------------------------------------------------------------
169  ! If RCI_REQUEST=0, then the solution was found with the required precision
170  !---------------------------------------------------------------------------
171  IF (rci_request.EQ.0) goto 3
172  !---------------------------------------------------------------------------
173  ! If RCI_REQUEST=1, then compute the vector A*TMP(IPAR(22))
174  ! and put the result in vector TMP(IPAR(23))
175  !---------------------------------------------------------------------------
176  IF (rci_request.EQ.1) THEN
177  CALL mkl_dcsrgemv('N',n, a, ia, ja, tmp(ipar(22)), tmp(ipar(23)))
178  goto 1
179  !---------------------------------------------------------------------------
180  ! If RCI_REQUEST=anything else, then DFGMRES subroutine failed
181  ! to compute the solution vector: COMPUTED_SOLUTION(N)
182  !---------------------------------------------------------------------------
183  ELSE
184  goto 999
185  ENDIF
186  !---------------------------------------------------------------------------
187  ! Reverse Communication ends here
188  ! Get the current iteration number and the FGMRES solution (DO NOT FORGET to
189  ! call DFGMRES_GET routine as COMPUTED_SOLUTION is still containing
190  ! the initial guess!)
191  !---------------------------------------------------------------------------
192 3 CALL dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)
193  !---------------------------------------------------------------------------
194  !---------------------------------------------------------------------------
195  ! Print solution vector: COMPUTED_SOLUTION(N) and
196  ! the number of iterations: ITERCOUNT
197  !---------------------------------------------------------------------------
198
199
200  !---------------------------------------------------------------------------
201  ! Release internal MKL memory that might be used for computations
202  ! NOTE: It is important to call the routine below to avoid memory leaks
203  ! unless you disable MKL Memory Manager
204  !---------------------------------------------------------------------------
205  CALL mkl_freebuffers
206
207  dvar=dnrm2(n,expected_solution,1)
208
209  !---------------------------------------------------------------------------
210  ! Release internal MKL memory that might be used for computations
211  ! NOTE: It is important to call the routine below to avoid memory leaks
212  ! unless you disable MKL Memory Manager
213  !---------------------------------------------------------------------------
214 999 WRITE( *,'(A,A,I5)') 'This example FAILED as the solver has',&
215  & ' returned the ERROR code', rci_request
216  CALL mkl_freebuffers
217  ! STOP 1
218  END SUBROUTINE fgmres_nopc
219
220  FUNCTION dnscsr(inmat,tolin)
221  REAL(sp),DIMENSION(:,:),INTENT(in)::inmat
222  REAL(sp),OPTIONAL::tolin
223  TYPE(sparse)::dnscsr
224  INTEGER::nz,m,n,job(6),info
225  REAL(sp)::tol
226  REAL(sp),DIMENSION(SIZE(inmat,1),SIZE(inmat,2))::mat
227  !! set the tolerance for zero search
228  IF(present(tolin)) THEN
229  tol=tolin
230  ELSE
231  tol=1.0e-8_sp
232  END IF
233  mat=inmat
234  !Zero out the entries that are below tolerance level
235  WHERE(abs(inmat)<tol)
236  mat=0.0_sp
237  END WHERE
238  ! Count the number of nonzero entries in the incoming matrix
239  nz=sparscalc(mat,tol)
240  !Number of rows in the incoming matrix
241  m=SIZE(mat,1)
242  !Number of columns in the incoming matrix
243  n=SIZE(mat,2)
244  ! Allocate the arrays for sparse storage
245  ALLOCATE(dnscsr%A(nz),dnscsr%IA(m+1),dnscsr%JA(nz))
246
247  ! Set the job preference
248  job(1)=0 ! Convert regular matrix to CSR
249  job(2)=1 ! Use one-based matrix notation
250  job(3)=1 ! Use one-based indexing in CSR format
251  job(4)=2 ! Input the whole matrix
252  job(5)=nz! Maximum number of nonzero elements
253  job(6)=1 ! Return Ia, Ja, and AA
254
255  CALL mkl_ddnscsr(job,m,n,mat,m,dnscsr%A,dnscsr%JA,dnscsr%IA,info)
256
257  END FUNCTION dnscsr
258
259  FUNCTION sparscalc(mat,tolin)
260  ! This function calculates the number of nonzero entries in
261  ! a sparse matrix. By default, numbers smaller than 1.0e-8
262  ! are considered zero, but the default can be overriden by
263  ! supplying the optional second paramter
264
265  IMPLICIT NONE
266  REAL(sp),DIMENSION(:,:),INTENT(in)::mat
267  REAL(sp),OPTIONAL::tolin
268  REAL(sp)::tol
269  INTEGER::nn,sparscalc
270  INTEGER,DIMENSION(SIZE(mat,1),SIZE(mat,2))::tmp
271  tmp=0
272  IF(present(tolin)) THEN
273  tol=tolin
274  ELSE
275  tol=1.0e-8_sp
276  END IF
277
278  WHERE(abs(mat)>tol)
279  tmp=1
280  END WHERE
281  sparscalc=sum(tmp)
282  END FUNCTION sparscalc
283
284  FUNCTION inverse(MAT)
285  IMPLICIT NONE
286  REAL(sp),DIMENSION(:,:),INTENT(in)::mat
287  REAL(SP),DIMENSION(SIZE(MAT,1),SIZE(MAT,2))::a,inverse
288  !This is a wrapper for inverting a matrix
289  INTEGER(I4B)::n,lda,info,lwork
290  INTEGER(I4B),DIMENSION(SIZE(MAT,1))::ipiv
291  REAL(sp),DIMENSION(:),ALLOCATABLE::work
292  n=SIZE(mat,1)
293  lda=n
294  inverse=mat
295  a=mat
296  !LU factorize the matrix first
297  CALL dgetrf(n,n,a,lda,ipiv,info)
298  lwork=n*64
299  ALLOCATE(work(lwork))
300  !Then invert
301  CALL dgetri( n, a, lda, ipiv, work, lwork, info )
302  DEALLOCATE(work)
303  inverse=a
304  print*,'inverse: maximum error',maxval(eye(n)-matmul(mat,a))
305
306  END FUNCTION inverse
307
308  SUBROUTINE eigen_val(A,wr,wi)
309  IMPLICIT NONE
310  REAL(sp),DIMENSION(:,:),INTENT(in)::a
311  REAL(sp),DIMENSION(SIZE(A,1)),INTENT(out)::wr,wi
312  REAL(sp),DIMENSION(SIZE(A,1))::scale,work,w
313  REAL(sp),DIMENSION(SIZE(A,1)-1)::tau
314  REAL(SP),DIMENSION(SIZE(A,1),SIZE(A,2))::z,h
315  INTEGER::ilo,ihi,info,n,lda,ldz,ldh,lwork
316  CHARACTER(len=1)::job,compz
317  job='B'
318  n=SIZE(a,1)
319  lda=n
320  lwork=n
321  ldz=n
322  ldh=n
323  h=a
324  ! Balance the incoming matrix
325  CALL dgebal(job, n, h, lda, ilo, ihi, scale, info)
326  ! Reduce the matrix to upper-Hysenberg form
327  CALL dgehrd(n, ilo, ihi, h, lda, tau, work, lwork, info)
328  ! Calculate eigenvalues
329  job='E' ! Only eigenvalues needed
330  compz='N' ! No Schur vectors are needed
331  z=0.0_sp ! It is not referenced when compz='N', but just set to zero
332  CALL dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
333  !hseqr(H, w, ilo,ihi,z,job,compz,info)
334  ! In the output w is the array of eigenvalues calculated by the routine
335
336  END SUBROUTINE eigen_val
337
338  FUNCTION eigentest(k)
339  IMPLICIT NONE
340  INTEGER,INTENT(in)::k
341  INTEGER::eigentest
342  REAL(sp),DIMENSION(8,8)::mat
343  REAL(sp),DIMENSION(8)::analytical,wr,wi,imag,rel
344
345
346  INTEGER::ii
347  mat(1,:)=(/611.0_sp,196.0_sp,-192.0_sp,407.0_sp,-8.0_sp,-52.0_sp,-49.0_sp,29.0_sp/)
348  mat(2,:)=(/196.0_sp,899.0_sp,113.0_sp,-192.0_sp,-71.0_sp,-43.0_sp,-8.0_sp,-44.0_sp/)
349  mat(3,:)=(/-192.0_sp,113.0_sp,899.0_sp,196.0_sp,61.0_sp,49.0_sp,8.0_sp,52.0_sp/)
350  mat(4,:)=(/407.0_sp,-192.0_sp,196.0_sp,611.0_sp,8.0_sp,44.0_sp,59.0_sp,-23.0_sp/)
351  mat(5,:)=(/-8.0_sp,-71.0_sp,61.0_sp,8.0_sp,411.0_sp,-599.0_sp,208.0_sp,208.0_sp/)
352  mat(6,:)=(/-52.0_sp,-43.0_sp,49.0_sp,44.0_sp,-599.0_sp,411.0_sp,208.0_sp,208.0_sp/)
353  mat(7,:)=(/-49.0_sp,-8.0_sp,8.0_sp,59.0_sp,208.0_sp,208.0_sp,99.0_sp,-911.0_sp/)
354  mat(8,:)=(/29.0_sp,-44.0_sp,52.0_sp,-23.0_sp,208.0_sp,208.0_sp,-911.0_sp,99.0_sp/)
355
356  print*,'eigentest: input matrix'
357  DO ii=1,8
358  WRITE(*,'(8(1x,f11.4))'),mat(ii,:)
359  END DO
360
361  CALL eigen_val(mat,wr,wi)
362
363
364  analytical=(/ 0.0_sp, 1020.0_sp, 510.0_sp-100.0_sp*sqrt(26.0_sp), &
365  & 510.0_sp+100.0_sp*sqrt(26.0_sp), 10.0_sp*sqrt(10405.0_sp), &
366  &-10.0_sp*sqrt(10405.0_sp), 1000.0_sp, 1000.0_sp/)
367
368
369  print*,'Calculated eigenvalues: analytical, real, imaginary'
370  DO ii=1,8
371  WRITE(*,'(3(1x,f11.4))'),analytical(ii),wr(ii),wi(ii)
372  END DO
373
374  mat(1,:)=(/10.0_sp,-65.0_sp,74.0_sp,42.0_sp,-22.0_sp,45.0_sp,51.0_sp,4.0_sp/)
375  mat(2,:)=(/25.0_sp,-69.0_sp,-18.0_sp,-18.0_sp,-9.0_sp,24.0_sp,-1.0_sp,79.0_sp/)
376  mat(3,:)=(/21.0_sp,-90.0_sp,45.0_sp,-64.0_sp,67.0_sp,-35.0_sp,49.0_sp,-42.0_sp/)
377  mat(4,:)=(/97.0_sp,80.0_sp,0.0_sp,74.0_sp,-77.0_sp,65.0_sp,86.0_sp,-72.0_sp/)
378  mat(5,:)=(/-79.0_sp,-97.0_sp,-65.0_sp,2.0_sp,0.0_sp,41.0_sp,47.0_sp,-5.0_sp/)
379  mat(6,:)=(/77.0_sp,-86.0_sp,69.0_sp,32.0_sp,-50.0_sp,0.0_sp,79.0_sp,95.0_sp/)
380  mat(7,:)=(/65.0_sp,-49.0_sp,-46.0_sp,-32.0_sp,-69.0_sp,88.0_sp,35.0_sp,-82.0_sp/)
381  mat(8,:)=(/3.0_sp,27.0_sp,-25.0_sp,3.0_sp,-14.0_sp,42.0_sp,-38.0_sp,-3.0_sp/)
382
383  print*,'eigentest: input matrix'
384  DO ii=1,8
385  WRITE(*,'(8(1x,f11.4))'),mat(ii,:)
386  END DO
387
388  CALL eigen_val(mat,wr,wi)
389
390  rel=(/140.01_sp,140.01_sp,-127.04_sp,-50.56_sp,&
391  &-50.56_sp,19.66_sp,19.66_sp,0.8_sp/)
392
393  imag=(/55.68_sp,-55.68_sp,0.0_sp,48.15_sp,&
394  &-48.15_sp,69.28_sp,-69.28_sp,0.0_sp/)
395
396  print*,'Eigenvalues: real(analytical), imaginary(analytical), real, imaginary'
397  DO ii=1,8
398  WRITE(*,'(4(1x,f11.4))'),rel(ii),imag(ii),wr(ii),wi(ii)
399  END DO
400
401
402  eigentest=1
403  END FUNCTION eigentest
404
405  FUNCTION filename(f,str)
406  INTEGER(I4B),INTENT(in)::f
407  CHARACTER(len=3),INTENT(in)::str
408  CHARACTER(len=40)::frmt
409  CHARACTER(len=18)::filename
410
411  IF(f<10) THEN
412  frmt='("../data/",a3,"00",i1,".vtk")'
413  WRITE(filename,frmt),str,f
414  elseif((f<100).AND.(f>=10)) THEN
415  frmt='("../data/",a3,"0",i2,".vtk")'
416  WRITE(filename,frmt),str,f
417  elseif((f<1000).AND.(f>=100)) THEN
418  frmt='("../data/",a3,i3,".vtk")'
419  WRITE(filename,frmt),str,f
420  END IF
421  END FUNCTION filename
422
423
424
425
426
427 END MODULE global