PROGRESS  master
prg_genz_mod.F90
Go to the documentation of this file.
1 
6 module prg_genz_mod
7 
8  use prg_openfiles_mod
9  use bml
10  use prg_kernelparser_mod
11  use prg_parallel_mod
12  use prg_extras_mod
13 
14  implicit none
15 
16  private !Everything is private by default
17 
18  integer, parameter :: dp = kind(1.0d0)
19 
20  public :: prg_buildzdiag, prg_buildzsparse, prg_parse_zsp, prg_init_zspmat
21  public :: prg_genz_sp_initial_zmat, prg_genz_sp_ref, prg_genz_sp_initialz0
22 
26  type, public :: genzspinp
27 
29  integer :: verbose
30 
32  integer :: nfirst
33 
35  integer :: nrefi
36 
38  integer :: nreff
39 
41  real(dp) :: numthresi
42 
44  real(dp) :: numthresf
45 
47  logical :: integration
48 
50  integer :: igenz
51 
53  logical :: zsp
54 
56  integer :: mdim
57 
59  character(20) :: bml_type
60 
61  end type genzspinp
62 
63 contains
64 
67  subroutine prg_parse_zsp(input,filename)
68  type(genzspinp), intent(inout) :: input
69  integer, parameter :: nkey_char = 1, nkey_int = 5, nkey_re = 2, nkey_log = 2
70  character(len=*) :: filename
71 
72  !Library of keywords with the respective defaults.
73  character(len=50), parameter :: keyvector_char(nkey_char) = [character(len=100) :: &
74  'BMLType=']
75  character(len=100) :: valvector_char(nkey_char) = [character(len=100) :: &
76  'Dense']
77 
78  character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: &
79  'Verbose=','NFirst=','NRefI=','NRefF=','MDim=']
80  integer :: valvector_int(nkey_int) = (/ &
81  0,10,3,1,-1 /)
82 
83  character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
84  'NumthreshI=','NumthreshF=' ]
85  real(dp) :: valvector_re(nkey_re) = (/&
86  1.0e-8, 1.0e-5 /)
87 
88  character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=100) :: &
89  'ZSP=','Int=']
90  logical :: valvector_log(nkey_log) = (/&
91  .false., .true. /)
92 
93  !Start and stop characters
94  character(len=50), parameter :: startstop(2) = [character(len=50) :: &
95  'ZSP{', '}']
96 
97  call prg_parsing_kernel(keyvector_char,valvector_char&
98  ,keyvector_int,valvector_int,keyvector_re,valvector_re,&
99  keyvector_log,valvector_log,trim(filename),startstop)
100 
101  !Characters
102  if(valvector_char(1) == "Dense")then
103  input%bml_type = bml_matrix_dense
104  elseif(valvector_char(1) == "Ellpack")then
105  input%bml_type = bml_matrix_ellpack
106  elseif(valvector_char(1) == "Ellblock")then
107  input%bml_type = bml_matrix_ellblock
108  endif
109 
110  !Integers
111  input%verbose = valvector_int(1)
112  input%nfirst = valvector_int(2)
113  input%nrefi = valvector_int(3)
114  input%nreff = valvector_int(4)
115  input%mdim = valvector_int(5)
116 
117  !Logicals
118  input%zsp = valvector_log(1)
119  input%integration = valvector_log(2)
120 
121  !Reals
122  input%numthresi = valvector_re(1)
123  input%numthresf = valvector_re(2)
124 
125  end subroutine prg_parse_zsp
126 
133  subroutine prg_init_zspmat(igenz,zk1_bml,zk2_bml,zk3_bml&
134  ,zk4_bml,zk5_bml,zk6_bml,norb,bml_type,bml_element_type)
135  integer :: norb, igenz
136  character(20) :: bml_type, my_bml_element_type
137  character(20), optional :: bml_element_type
138  type(bml_matrix_t) :: zk1_bml
139  type(bml_matrix_t) :: zk2_bml
140  type(bml_matrix_t) :: zk3_bml
141  type(bml_matrix_t) :: zk4_bml
142  type(bml_matrix_t) :: zk5_bml
143  type(bml_matrix_t) :: zk6_bml
144 
145  if(present(bml_element_type))then
146  my_bml_element_type = bml_element_type
147  else
148  my_bml_element_type = bml_element_real
149  endif
150 
151  igenz = 0
152 
153  if(bml_get_n(zk1_bml).le.0)then
154  call bml_zero_matrix(bml_type,my_bml_element_type,dp,norb,norb,zk1_bml)
155  call bml_zero_matrix(bml_type,my_bml_element_type,dp,norb,norb,zk2_bml)
156  call bml_zero_matrix(bml_type,my_bml_element_type,dp,norb,norb,zk3_bml)
157  call bml_zero_matrix(bml_type,my_bml_element_type,dp,norb,norb,zk4_bml)
158  call bml_zero_matrix(bml_type,my_bml_element_type,dp,norb,norb,zk5_bml)
159  call bml_zero_matrix(bml_type,my_bml_element_type,dp,norb,norb,zk6_bml)
160  endif
161 
162  end subroutine prg_init_zspmat
163 
174  subroutine prg_buildzdiag(smat_bml,zmat_bml,threshold,mdimin,bml_type,verbose)
175  ! use extras
176  real(dp) :: err_check
177  character(len=*), intent(in) :: bml_type
178  character(20) :: bml_element_type
179  integer :: i, j, mdim, norb
180  integer, intent(in) :: mdimin
181  integer, optional, intent(in) :: verbose
182  real(dp) :: mls_i
183  real(dp) :: invsqrt, threshold
184  real(dp), allocatable :: nono_evals(:), nonotmp(:,:), smat(:,:), umat(:,:)
185  real(dp), allocatable :: zmat(:,:)
186  type(bml_matrix_t) :: nonotmp_bml, saux_bml, umat_bml, umat_t_bml
187  type(bml_matrix_t) :: zmat_bml
188  type(bml_matrix_t), intent(inout) :: smat_bml
189 
190  if(present(verbose))then
191  if(verbose >= 1)then
192  write(*,*)""; write(*,*)"In buildzdiag ..."
193  endif
194  endif
195 
196  if(bml_get_precision(smat_bml) == 1 .or.&
197  &bml_get_precision(smat_bml) == 2)then
198  bml_element_type = "real"
199  elseif(bml_get_precision(smat_bml) == 3 .or.&
200  &bml_get_precision(smat_bml) == 4)then
201  bml_element_type = "complex"
202  endif
203 
204  norb = bml_get_n(smat_bml)
205  mdim = bml_get_m(smat_bml)
206 
207  !Allocate temporary matrices.
208  allocate(nono_evals(norb))
209  allocate(umat(norb,norb))
210  allocate(nonotmp(norb,norb))
211 
212  if(.not. bml_allocated(zmat_bml))then
213  call bml_zero_matrix(bml_matrix_dense,bml_element_type,dp,norb,norb,zmat_bml)
214  endif
215 
216  !Auxiliary matrices.
217  call bml_zero_matrix(bml_type,bml_element_type,dp,norb,norb,umat_bml)
218  call bml_zero_matrix(bml_type,bml_element_type,dp,norb,norb,nonotmp_bml)
219 
220  !Eigenvectors and eigenalues of the overlap s.
221  mls_i = mls()
222  call bml_diagonalize(smat_bml,nono_evals,umat_bml)
223 
224  if(present(verbose))then
225  if(verbose >= 1)then
226  write(*,*)"Time for S diag = "//to_string(mls() - mls_i)//" ms"
227  endif
228  endif
229  call bml_export_to_dense(umat_bml, umat)
230 
231  nonotmp=0.0_dp
232 
233  !Doing s^-1/2 u^t
234 
235  !$omp parallel do default(none) &
236  !$omp shared(norb,nono_evals,umat,nonotmp) &
237  !$omp private(i,j,invsqrt)
238  do i = 1, norb
239  if(nono_evals(i).lt.0.0_dp) stop 'matrix s has a 0 eigenvalue'
240  invsqrt = 1.0_dp/sqrt(nono_evals(i))
241  do j = 1, norb
242  nonotmp(i,j) = umat(j,i) * invsqrt
243  end do
244  end do
245  !$omp end parallel do
246 
247  call bml_import_from_dense(bml_type, nonotmp, nonotmp_bml, threshold, mdim)
248 
249  !Doing u s^-1/2 u^t
250  call bml_multiply(umat_bml,nonotmp_bml,zmat_bml,1.0_dp,0.0_dp,threshold)
251 
252  deallocate(nonotmp)
253  deallocate(nono_evals)
254  deallocate(umat)
255  call bml_deallocate(umat_bml)
256  call bml_deallocate(nonotmp_bml)
257 
258  end subroutine prg_buildzdiag
259 
272  subroutine prg_buildzsparse(smat_bml,zmat_bml,igenz,mdim,&
273  bml_type,zk1_bml,zk2_bml,zk3_bml&
274  ,zk4_bml,zk5_bml,zk6_bml,nfirst,nrefi,nreff&
275  ,thresholdi,thresholdf,integration,verbose)
276  !use extras
277  !real(dp) :: err_check
278  !real(dp), allocatable :: smat(:,:), zmat(:,:)
279  character(20) :: bml_type
280  integer :: kk, i, igenz, mdim
281  integer :: nfirst, norb, nref, nreff
282  integer :: nrefi, verbose
283  logical :: integration
284  real(dp) :: sec_i, sec_ii
285  real(dp) :: threshold, thresholdf, thresholdi
286  type(bml_matrix_t) :: smat_bml, zk1_bml, zk2_bml, zk3_bml
287  type(bml_matrix_t) :: zk4_bml, zk5_bml, zk6_bml, zmat_bml
288 
289  norb = bml_get_n(smat_bml)
290  if(mdim<0) mdim = norb
291 
292  kk=6 !Total number of stored z matrices for the xl integration scheme.
293 
294  if(verbose.eq.1)sec_ii=mls() !Gets the actual time in ms.
295 
296  ! Number of refinements and threshold value for the first "nfirst" iterations.
297  if(igenz.lt.nfirst)then
298  nref=nrefi; threshold = thresholdi !For the firsts iterations.
299  else
300  nref=nreff; threshold = thresholdf !For the following iterations.
301  end if
302 
303  if(verbose.eq.1)sec_i=mls() !Firs calculation of z using the graph approach.
304  if(igenz.eq.1) call prg_genz_sp_initialz0(smat_bml,zmat_bml,norb,mdim,bml_type,threshold)
305  if(verbose.eq.1)write(*,*)"Time for prg_initial estimate "//to_string(mls()-sec_i)//" ms"
306 
307  if(verbose.eq.1)sec_i=mls()! integration scheme.
308  if(integration)then
309  call prg_genz_sp_int(zmat_bml,zk1_bml,zk2_bml,zk3_bml&
310  &,zk4_bml,zk5_bml,zk6_bml,igenz,norb,bml_type,threshold)
311  end if
312  if(verbose.eq.1)write(*,*)"Time for xl scheme "//to_string(mls()-sec_i)//" ms"
313 
314  if(verbose.eq.1)sec_i=mls()! Refinement.
315  call prg_genz_sp_ref(smat_bml,zmat_bml,nref,mdim,bml_type,threshold)
316  if(verbose.eq.1)write(*,*)"Time for prg_genz_sp_ref "//to_string(mls()-sec_i)//" ms"
317 
318  ! call bml_export_to_dense()
319  ! call prg_delta(xmat,smat,norb,err_check) !to check for the accuracy of the approximation (prg_delta)
320  ! write(*,*)"err", err_check, norb
321  ! stop
322 
323  if(verbose.eq.1)write(*,*)"Time for prg_buildZsparse "//to_string(igenz)//" "//to_string(mls()-sec_ii)//" ms"
324 
325  end subroutine prg_buildzsparse
326 
338  subroutine prg_genz_sp_initialz0(smat_bml,zmat_bml,norb,mdim,bml_type_f,threshold)
339  character(20) :: bml_type, bml_type_f, bml_element_type
340  integer :: i, ii, j, jj
341  integer :: k, l, mdim, norb
342  integer, allocatable :: kk(:)
343  real(dp) :: err_check, invsqrt, threshold, thresholdz0
344  real(dp), allocatable :: sitmp(:,:), smat(:,:), sthres(:,:), stmp(:,:)
345  real(dp), allocatable :: stmp_evals(:), utmp(:,:), zmat(:,:), ztmp(:,:)
346  type(bml_matrix_t) :: sitmp_bml, sthres_bml, stmp_bml, utmp_bml
347  type(bml_matrix_t) :: xtmp_bml, zmat_bml
348  type(bml_matrix_t), intent(in) :: smat_bml
349 
350  !When this subroutine is parallelized, the following threshold can be set to larger (coarser) values.
351  !The reason for this is to reduce the memory per thread.
352  thresholdz0= 1.0d-10 !This could be passed as an input.
353 
354  bml_type= bml_matrix_dense !All the operations are performed in bml_dense.
355  ! bml_element_type = bml_get_element_type(smat_bml)
356  if(bml_get_precision(smat_bml) == 1 .or.&
357  &bml_get_precision(smat_bml) == 2)then
358  bml_element_type = "real"
359  elseif(bml_get_precision(smat_bml) == 3 .or.&
360  &bml_get_precision(smat_bml) == 4)then
361  bml_element_type = "complex"
362  endif
363 
364  !Part of the operations are still done in pure dense format.
365  allocate(zmat(norb,norb))
366  zmat = 0.0_dp
367 
368  allocate(sthres(norb,norb))
369 
370  allocate(smat(norb,norb))
371  call bml_export_to_dense(smat_bml,smat)
372 
373  call bml_zero_matrix(bml_type,bml_element_type,dp,norb,norb,sthres_bml) !a thresholded version of s
374 
375  call bml_import_from_dense(bml_type,smat,sthres_bml,thresholdz0,mdim)
376 
377  call bml_export_to_dense(sthres_bml,sthres) !We will use the dense to extract the small blocks.
378 
379  allocate(kk(norb)) !The nonzero positions in the row
380 
381  !NOTE: The following loop must be OMP parallelized as parallelization is trivial, however, as BML
382  !operations are taking all the threads this is not possible. BML needs to recognize automatically
383  !when "threads" are requested by the hosting code.
384  do i = 1, norb !Z0 is the prg_initial guess for the iterative refinement.
385 
386  jj=0
387  kk=0
388  do j=1,norb
389  if(abs(sthres(i,j)).gt.thresholdz0)then
390  jj=jj+1
391  kk(jj)=j
392  endif
393  enddo
394  k=jj
395 
396  !The followings will be dense matrices since they are the small blocks.
397  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,stmp_bml)
398  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,sitmp_bml)
399  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,utmp_bml)
400  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,xtmp_bml)
401 
402  allocate(stmp(k,k),sitmp(k,k))
403  allocate(utmp(k,k),stmp_evals(k),ztmp(k,k))
404 
405  stmp = 0.0_dp !Small dense block to be extracted from s.
406  stmp_evals = 0.0_dp !eigenvalues for the small dense s matrices.
407  sitmp = 0.0_dp
408  utmp = 0.0_dp !Matrix containing the eigenvectors of small s
409  ztmp = 0.0_dp
410 
411  do j = 1, k !Extracting the small s matrices
412  do l = 1, k
413  stmp(l,j) = smat(kk(l),kk(j))
414  end do
415  end do
416 
417  call bml_import_from_dense(bml_type, stmp, stmp_bml)
418  call bml_diagonalize(stmp_bml,stmp_evals,utmp_bml)
419  call bml_export_to_dense(utmp_bml,utmp)
420 
421  do j = 1, k !Applying the inverse sqrt function to the eigenvalues.
422  invsqrt = 1.0_dp/sqrt(stmp_evals(j))
423  do l = 1, k
424  sitmp(l,j) = utmp(l,j) * invsqrt
425  end do
426  end do
427 
428  utmp=transpose(utmp)
429 
430  call bml_import_from_dense(bml_type , utmp , utmp_bml)
431  call bml_import_from_dense(bml_type, sitmp, sitmp_bml)
432  call bml_multiply(sitmp_bml,utmp_bml,xtmp_bml,1.0_dp,0.0_dp)
433  call bml_export_to_dense(xtmp_bml,ztmp)
434 
435  do l = 1, k !Reconstructing the large Z0 based in the small Z.
436  do j = 1, k
437  if(kk(j).eq.i)then !For the row corresponding to the extracted dense block.
438  zmat(i,kk(l)) = ztmp(j,l)
439  end if
440  end do
441  end do
442 
443  deallocate(stmp,sitmp) !Deallocate temporary matrices
444  deallocate(utmp)
445  deallocate(stmp_evals)
446  deallocate(ztmp)
447  call bml_deallocate(sitmp_bml)
448  call bml_deallocate(stmp_bml)
449  call bml_deallocate(utmp_bml)
450 
451  end do
452 
453  call bml_zero_matrix(bml_type_f,bml_element_type,dp,norb,norb,zmat_bml)
454  call bml_import_from_dense(bml_type_f,zmat, zmat_bml,threshold,mdim) !Converting to bml format
455 
456  ! call prg_delta(zmat,smat,norb,err_check) !to check for the accuracy of the approximation (prg_delta)
457  ! call sparsity(smat,norb,spa)
458  ! write(*,*)"err", err_check, norb
459  ! stop
460 
461  deallocate(sthres)
462  deallocate(kk)
463  call bml_deallocate(sthres_bml)
464 
465  end subroutine prg_genz_sp_initialz0
466 
478  subroutine prg_genz_sp_initial_zmat(smat_bml,zmat_bml,norb,mdim,bml_type_f,threshold)
479  ! use extras
480  implicit none
481  character(20) :: bml_type, bml_type_f, bml_element_type
482  integer :: i, ii, j, jj
483  integer :: k, l, mdim, norb
484  integer, allocatable :: kk(:)
485  real(dp), intent(in) :: threshold
486  real(dp) :: err_check, invsqrt
487  real(dp) :: thresholdz0
488  real(dp), allocatable :: sitmp(:,:), smat(:,:), sthres(:,:)
489  real(dp), allocatable :: stmp(:,:)
490  real(dp), allocatable :: stmp_evals(:), utmp(:,:)
491  real(dp), allocatable :: zmat(:,:), ztmp(:,:)
492  type(bml_matrix_t) :: sitmp_bml, sthres_bml, stmp_bml
493  type(bml_matrix_t) :: xtmp_bml, zmat_bml, utmp_bml
494  type(bml_matrix_t), intent(in) :: smat_bml
495 
496  !When this subroutine is parallelized, the following threshold can be set to
497  !larger (coarser) values.
498  !The reason for this is to reduce the memory per thread.
499  !thresholdz0= 1.0d-10 !This could be passed as an input.
500 
501  bml_type= bml_matrix_dense !All the operations are performed in bml_dense.
502  ! bml_element_type = bml_get_element_type(smat_bml)
503  if(bml_get_precision(smat_bml) == 1 .or.&
504  &bml_get_precision(smat_bml) == 2)then
505  bml_element_type = "real"
506  elseif(bml_get_precision(smat_bml) == 3 .or.&
507  &bml_get_precision(smat_bml) == 4)then
508  bml_element_type = "complex"
509  endif
510 
511  !Part of the operations are still done in pure dense format.
512  allocate(zmat(norb,norb))
513  allocate(sthres(norb,norb))
514 
515  allocate(smat(norb,norb))
516  call bml_export_to_dense(smat_bml,smat)
517 
518  ! a thresholded version of s
519  call bml_zero_matrix(bml_type,bml_element_type,dp,norb,norb,sthres_bml)
520 
521  !call bml_import_from_dense(bml_type,smat,sthres_bml,thresholdZ0,mdim)
522  call bml_import_from_dense(bml_type,smat,sthres_bml,threshold,mdim)
523 
524  ! We will use the dense to extract the small blocks.
525  call bml_export_to_dense(sthres_bml,sthres)
526 
527  allocate(kk(norb)) !The nonzero positions in the row
528 
529  !NOTE: The following loop must be OMP parallelized as parallelization is
530  !trivial, however, as BML
531  !operations are taking all the threads this is not possible. BML needs to
532  !recognize automatically
533  !when "threads" are requested by the hosting code.
534  do i = 1, norb !Z0 is the prg_initial guess for the iterative refinement.
535  !do i = bml_getLocalRowMin(smat_bml, getMyRank()), &
536  ! bml_getLocalRowMax(smat_bml, getMyRank())
537 
538  jj=0
539  kk=0
540  do j=1,norb
541  !if(abs(sthres(i,j)).gt.thresholdZ0)then
542  if(abs(sthres(i,j)).gt.threshold)then
543  jj=jj+1
544  kk(jj)=j
545  endif
546  enddo
547  k=jj
548 
549  !The followings will be dense matrices since they are the small blocks.
550  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,stmp_bml)
551  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,sitmp_bml)
552  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,utmp_bml)
553  call bml_zero_matrix(bml_type,bml_element_type,dp, k,k,xtmp_bml)
554 
555  allocate(stmp(k,k),sitmp(k,k))
556  allocate(utmp(k,k),stmp_evals(k),ztmp(k,k))
557 
558  stmp = 0.0_dp !Small dense block to be extracted from s.
559  stmp_evals = 0.0_dp !eigenvalues for the small dense s matrices.
560  sitmp = 0.0_dp
561  utmp = 0.0_dp !Matrix containing the eigenvectors of small s
562  ztmp = 0.0_dp
563 
564  do j = 1, k !Extracting the small s matrices
565  do l = 1, k
566  stmp(l,j) = smat(kk(l),kk(j))
567  end do
568  end do
569 
570  call bml_import_from_dense(bml_matrix_dense, stmp, stmp_bml)
571  call bml_diagonalize(stmp_bml,stmp_evals,utmp_bml)
572  call bml_export_to_dense(utmp_bml,utmp)
573 
574  do j = 1, k !Applying the inverse sqrt function to the eigenvalues.
575  invsqrt = 1.0_dp/sqrt(stmp_evals(j))
576  do l = 1, k
577  sitmp(l,j) = utmp(l,j) * invsqrt
578  end do
579  end do
580 
581  utmp=transpose(utmp)
582 
583  call bml_import_from_dense(bml_type, utmp, utmp_bml)
584  call bml_import_from_dense(bml_type, sitmp, sitmp_bml)
585  call bml_multiply(sitmp_bml,utmp_bml,xtmp_bml,1.0_dp,0.0_dp)
586  call bml_export_to_dense(xtmp_bml,ztmp)
587 
588  do l = 1, k !Reconstructing the large Z0 based in the small Z.
589  do j = 1, k
590  if(kk(j).eq.i)then !For the row corresponding to the extracted dense
591  !block.
592  zmat(i,kk(l)) = ztmp(j,l)
593  end if
594  end do
595  end do
596 
597  deallocate(stmp,sitmp) !Deallocate temporary matrices
598  deallocate(utmp)
599  deallocate(stmp_evals)
600  deallocate(ztmp)
601  call bml_deallocate(sitmp_bml)
602  call bml_deallocate(stmp_bml)
603  call bml_deallocate(utmp_bml)
604 
605  end do
606 
607  !call bml_zero_matrix(bml_type_f,bml_element_type,bml_element_precision,norb,norb,zmat_bml)
608  call bml_import_from_dense(bml_type_f,zmat,zmat_bml,threshold,mdim, &
609  bml_get_distribution_mode(smat_bml))
610  !Converting to bml format
611 
612  ! call prg_delta(zmat,smat,norb,err_check) !to check for the accuracy of
613  ! the approximation (prg_delta)
614  ! call sparsity(smat,norb,spa)
615  ! write(*,*)"err", err_check, norb
616  ! stop
617  deallocate(sthres)
618  deallocate(kk)
619  call bml_deallocate(sthres_bml)
620 
621  end subroutine prg_genz_sp_initial_zmat
622 
633  subroutine prg_genz_sp_int(zmat_bml,zk1_bml,zk2_bml,zk3_bml&
634  &,zk4_bml,zk5_bml,zk6_bml,igenz,norb,bml_type&
635  &,threshold)
636  integer :: igenz,norb,KK
637  character(20) :: bml_element_type
638  real(dp) :: alpha, kappa, c0, c1, c2, c3, c4, c5
639  real(dp) :: threshold
640  type(bml_matrix_t) :: zmat_bml
641  type(bml_matrix_t) :: zk1_bml,zk2_bml,zk3_bml,zk4_bml,zk5_bml,zk6_bml
642  character(20) :: bml_type
643 
644  kk=6
645  alpha=0.0180_dp
646  kappa=1.82_dp
647 
648  ! bml_element_type = bml_get_element_type(zmat_bml)
649  if(bml_get_precision(zmat_bml) == 1 .or.&
650  &bml_get_precision(zmat_bml) == 2)then
651  bml_element_type = "real"
652  elseif(bml_get_precision(zmat_bml) == 3 .or.&
653  &bml_get_precision(zmat_bml) == 4)then
654  bml_element_type = "complex"
655  endif
656 
657  !The following constants are the original constants premultiplied by alpha.
658 
659  c0=-0.1080_dp
660  c1=0.2520_dp
661  c2=-0.1440_dp
662  c3=-0.0540_dp
663  c4=0.0720_dp
664  c5=-0.0180_dp
665 
666 
667  if(igenz.eq.kk)then
668 
669  call bml_zero_matrix(bml_type,bml_element_type,dp,norb ,norb,zk1_bml)
670  call bml_zero_matrix(bml_type,bml_element_type,dp,norb ,norb,zk2_bml)
671  call bml_zero_matrix(bml_type,bml_element_type,dp,norb ,norb,zk3_bml)
672  call bml_zero_matrix(bml_type,bml_element_type,dp,norb ,norb,zk4_bml)
673  call bml_zero_matrix(bml_type,bml_element_type,dp,norb ,norb,zk5_bml)
674  call bml_zero_matrix(bml_type,bml_element_type,dp,norb ,norb,zk6_bml)
675 
676  call bml_copy(zmat_bml,zk1_bml);call bml_copy(zmat_bml,zk2_bml);call bml_copy(zmat_bml,zk3_bml)
677  call bml_copy(zmat_bml,zk4_bml);call bml_copy(zmat_bml,zk5_bml);call bml_copy(zmat_bml,zk6_bml)
678 
679  end if
680 
681  if(igenz.ge.kk)then !Here we change z by applying the t-r integration scheme
682 
683  call bml_add_deprecated(1.0_dp,zmat_bml,-1.0_dp,zk6_bml,threshold) !Z(t)-Z~(t)
684  call bml_scale(kappa,zmat_bml)
685  call bml_add_deprecated(1.0_dp,zmat_bml,2.0_dp,zk6_bml,threshold) !Z(t)+2Z~(t)
686  call bml_add_deprecated(1.0_dp,zmat_bml,-1.0_dp,zk5_bml,threshold) !Z(t)-Z~(t-dt)
687 
688  !Dissipation force term:
689  call bml_add_deprecated(1.0_dp,zmat_bml,c0,zk6_bml,threshold) !Z(t)+c0*Z~(t)
690  call bml_add_deprecated(1.0_dp,zmat_bml,c1,zk5_bml,threshold)
691  call bml_add_deprecated(1.0_dp,zmat_bml,c2,zk4_bml,threshold)
692  call bml_add_deprecated(1.0_dp,zmat_bml,c3,zk3_bml,threshold)
693  call bml_add_deprecated(1.0_dp,zmat_bml,c4,zk2_bml,threshold)
694  call bml_add_deprecated(1.0_dp,zmat_bml,c5,zk1_bml,threshold) !Z(t)+c0*Z~(t-5*dt)
695 
696  end if
697 
698  if(igenz.ge.kk)then !Here we are shifting the z matrices.
699 
700  call bml_copy(zk2_bml,zk1_bml) !Z~(t-5*dt)=Z~(t-4*dt)
701  call bml_copy(zk3_bml,zk2_bml) !Z~(t-4*dt)=Z~(t-3*dt)
702  call bml_copy(zk4_bml,zk3_bml)
703  call bml_copy(zk5_bml,zk4_bml)
704  call bml_copy(zk6_bml,zk5_bml) !Z~(t-dt)=Z~(t)
705  call bml_copy(zmat_bml,zk6_bml) !Z~(t)=Z~(t+dt)
706 
707  end if
708 
709  end subroutine prg_genz_sp_int
710 
719  subroutine prg_genz_sp_ref(smat_bml,zmat_bml,nref,norb,bml_type,threshold)
720  integer :: k
721  integer, intent(inout) :: norb
722  integer, intent(in) :: nref
723  real(dp) :: err_check, sec_i
724  real(dp), allocatable :: smat(:,:), zmat(:,:)
725  real(dp),intent(in) :: threshold
726  type(bml_matrix_t) :: temp_bml
727  type(bml_matrix_t) :: temp1_bml
728  type(bml_matrix_t) :: temp2_bml
729  type(bml_matrix_t) :: idscaled_bml
730  type(bml_matrix_t) :: xmat_t_bml
731  type(bml_matrix_t),intent(inout) :: zmat_bml
732  type(bml_matrix_t) :: aux_bml
733  type(bml_matrix_t), intent(in) :: smat_bml
734  character(20),intent(in) :: bml_type
735  character(20) :: bml_element_type
736 
737  norb = bml_get_n(smat_bml)
738 
739  ! bml_element_type = bml_get_element_type(smat_bml)
740  if(bml_get_precision(smat_bml) == 1 .or.&
741  &bml_get_precision(smat_bml) == 2)then
742  bml_element_type = "real"
743  elseif(bml_get_precision(smat_bml) == 3 .or.&
744  &bml_get_precision(smat_bml) == 4)then
745  bml_element_type = "complex"
746  endif
747 
748  call bml_zero_matrix(bml_type,bml_element_type,dp,norb,norb,idscaled_bml)
749 
750  call bml_add_identity(idscaled_bml, 1.0_dp, threshold) !1.0 [0] + 1.0 I
751  call bml_scale(1.8750_dp,idscaled_bml) ! 1.875*I
752 
753  call bml_noinit_matrix(bml_type,bml_element_type,dp,norb ,norb,temp_bml)
754  call bml_noinit_matrix(bml_type,bml_element_type,dp,norb ,norb,temp2_bml)
755 
756  sec_i=mls() !Firs calculation of z using the graph approach.
757  do k = 1, nref !Iterative refinement
758 
759  !Enforcing symmetry (in bml).
760  call bml_transpose(zmat_bml, xmat_t_bml) !Z^t
761  call bml_add_deprecated(0.50_dp,zmat_bml, 0.50_dp, xmat_t_bml,threshold) !(Z^t+Z)/2
762 
763  call bml_multiply(smat_bml,zmat_bml,temp_bml, 1.0_dp, 0.0_dp,threshold) !S*Z
764 
765  call bml_multiply(zmat_bml,temp_bml, temp2_bml, 1.0_dp, 0.0_dp,threshold) !X = Z^t*S*Z
766 
767  call bml_multiply(temp2_bml, temp2_bml, temp_bml, 1.0_dp, 0.0_dp,threshold) !X*X
768 
769  call bml_scale(0.3750_dp, temp_bml)
770  call bml_scale(-1.250_dp, temp2_bml)
771 
772  !Temp = 1.875*I - 1.25*X + 0.375*X^2
773  call bml_add_deprecated(1.0_dp,temp2_bml,1.0_dp, idscaled_bml,threshold)
774  call bml_add_deprecated(1.0_dp,temp_bml,1.0_dp, temp2_bml,threshold)
775 
776  call bml_multiply(zmat_bml,temp_bml,temp2_bml, 1.0_dp, 0.0_dp,threshold) !Z*Temp
777 
778  call bml_copy(temp2_bml,zmat_bml)
779 
780  end do
781  call bml_deallocate(temp_bml)
782  call bml_deallocate(temp1_bml)
783  call bml_deallocate(temp2_bml)
784 
785  write(*,*)"Time for ref loop "//to_string(mls()-sec_i)//" ms"
786 
787 
788  call bml_deallocate(idscaled_bml)
789  call bml_deallocate(xmat_t_bml)
790 
791  end subroutine prg_genz_sp_ref
792 
793 end module prg_genz_mod
prg_extras_mod::to_string
Definition: prg_extras_mod.F90:25
prg_extras_mod
Extra routines.
Definition: prg_extras_mod.F90:6
prg_extras_mod::mls
real(dp) function, public mls()
To get the actual time in milliseconds.
Definition: prg_extras_mod.F90:140
prg_extras_mod::dp
integer, parameter dp
Definition: prg_extras_mod.F90:31
prg_genz_mod
To produce a matrix which is needed to orthogonalize .
Definition: prg_genz_mod.F90:6
prg_genz_mod::prg_buildzsparse
subroutine, public prg_buildzsparse(smat_bml, zmat_bml, igenz, mdim, bml_type, zk1_bml, zk2_bml, zk3_bml, zk4_bml, zk5_bml, zk6_bml, nfirst, nrefi, nreff, thresholdi, thresholdf, integration, verbose)
Inverse factorization using Niklasson's algorithm.
Definition: prg_genz_mod.F90:276
prg_genz_mod::prg_genz_sp_initial_zmat
subroutine, public prg_genz_sp_initial_zmat(smat_bml, zmat_bml, norb, mdim, bml_type_f, threshold)
Initial estimation of Z.
Definition: prg_genz_mod.F90:479
prg_genz_mod::prg_buildzdiag
subroutine, public prg_buildzdiag(smat_bml, zmat_bml, threshold, mdimin, bml_type, verbose)
Usual subroutine involving diagonalization. , where = eigenvectors and = eigenvalues....
Definition: prg_genz_mod.F90:175
prg_genz_mod::prg_genz_sp_ref
subroutine, public prg_genz_sp_ref(smat_bml, zmat_bml, nref, norb, bml_type, threshold)
Iterative refinement.
Definition: prg_genz_mod.F90:720
prg_genz_mod::prg_parse_zsp
subroutine, public prg_parse_zsp(input, filename)
The parser for genz solver.
Definition: prg_genz_mod.F90:68
prg_genz_mod::prg_genz_sp_int
subroutine prg_genz_sp_int(zmat_bml, zk1_bml, zk2_bml, zk3_bml, zk4_bml, zk5_bml, zk6_bml, igenz, norb, bml_type, threshold)
Inverse factorization using Niklasson's algorithm.
Definition: prg_genz_mod.F90:636
prg_genz_mod::prg_genz_sp_initialz0
subroutine, public prg_genz_sp_initialz0(smat_bml, zmat_bml, norb, mdim, bml_type_f, threshold)
Initial estimation of Z.
Definition: prg_genz_mod.F90:339
prg_genz_mod::prg_init_zspmat
subroutine, public prg_init_zspmat(igenz, zk1_bml, zk2_bml, zk3_bml, zk4_bml, zk5_bml, zk6_bml, norb, bml_type, bml_element_type)
Initiates the matrices for the Xl integration of Z.
Definition: prg_genz_mod.F90:135
prg_kernelparser_mod
Some general parsing functions.
Definition: prg_kernelparser_mod.F90:4
prg_kernelparser_mod::prg_parsing_kernel
subroutine, public prg_parsing_kernel(keyvector_char, valvector_char, keyvector_int, valvector_int, keyvector_re, valvector_re, keyvector_log, valvector_log, filename, startstop)
The general parsing function. It is used to vectorize a set of "keywords" "value" pairs as included i...
Definition: prg_kernelparser_mod.F90:30
prg_openfiles_mod
Module to handle input output files for the PROGRESS lib.
Definition: prg_openfiles_mod.F90:4
prg_parallel_mod
The parallel module.
Definition: prg_parallel_mod.F90:6
prg_genz_mod::genzspinp
Input for the genz driver. This type controlls all the variables that are needed by genz.
Definition: prg_genz_mod.F90:26