Line data Source code
1 : !!****m* ABINIT/m_exp_mat
2 : !! NAME
3 : !! m_exp_mat
4 : !!
5 : !! FUNCTION
6 : !! This subroutine calculate the exponential of a matrix
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2002-2026 ABINIT group (XG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
14 : !!
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : MODULE m_exp_mat
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_errors
29 : use m_linalg_interfaces
30 :
31 : implicit none
32 :
33 : private
34 :
35 : public :: exp_mat ! exponential of a complex matrix
36 :
37 :
38 : interface exp_mat
39 : module procedure exp_mat_cx
40 : end interface exp_mat
41 :
42 :
43 : CONTAINS !===========================================================
44 : !!***
45 :
46 : !!****f* m_exp_mat/exp_mat_cx
47 : !! NAME
48 : !! exp_mat_cx
49 : !!
50 : !! FUNCTION
51 : !! Returns the exponential of a complex matrix multiplied a scalar
52 : !!
53 : !! INPUTS
54 : !! mat_a = the complex matrix
55 : !! mat_a = the size of mat_a
56 : !! factor = a real factor multiplying at the exponent
57 : !!
58 : !! OUTPUT
59 : !! exp_mat_cx = its exponential is returned in the same matrix
60 : !! SOURCE
61 :
62 0 : subroutine exp_mat_cx(mat_a,mat_a_size,factor)
63 :
64 : !Arguments ------------------------------------
65 : ! scalars
66 : real(dp),intent(in) :: factor
67 : ! arrays
68 : complex(dp),intent(inout) :: mat_a(:,:)
69 : !complex(dp) :: exp_mat_cx(mat_a_size,mat_a_size)
70 :
71 : !Local ------------------------------------------
72 : ! scalars
73 : integer :: info,mat_a_size,ii
74 : integer,parameter :: maxsize=3
75 : integer,parameter :: lwork=(1+32)*maxsize
76 :
77 : ! arrays
78 0 : integer :: ipvt(mat_a_size)
79 : character(len=500) :: msg
80 : real(dp) :: rwork(2*maxsize)
81 0 : complex(dp),allocatable :: ww(:),uu(:,:)
82 : complex(dp) :: work(lwork),vl(1,1)
83 : ! *********************************************************************
84 :
85 0 : mat_a_size = max(1,size(mat_a,dim=1))
86 0 : ABI_MALLOC(ww,(mat_a_size))
87 0 : ABI_MALLOC(uu,(mat_a_size,mat_a_size))
88 :
89 : !Now it calculates the eigenvalues and eigenvectors of the matrix
90 : call ZGEEV('No left vectors','Vectors (right)',mat_a_size, mat_a, mat_a_size,ww,&
91 0 : vl,1,uu, mat_a_size, work, lwork, rwork, info)
92 0 : if (info/=0) then
93 0 : write(msg,'(a,i4)')'Wrong value for rwork ',info
94 0 : ABI_BUG(msg)
95 : end if
96 :
97 : !!debbug
98 : ! write(std_out,*)'mat_a',mat_a
99 : ! write(std_out,*)'mat_a_size',mat_a_size
100 :
101 : ! write(std_out,*)'eigenvalues'
102 : ! write(std_out,*)ww
103 : ! write(std_out,*)'eigenvectors'
104 : ! do ii=1,g_mat_size
105 : ! write(std_out,*)'autov',ii
106 : ! write(std_out,*)uu(:,ii)
107 : ! end do
108 : ! write(std_out,*)'optimal workl=',work(1)
109 :
110 : ! write(std_out,*)'info', info
111 : ! write(std_out,*) '------------------'
112 : !!enddebug
113 :
114 :
115 :
116 : !-----------------------------------------------------------
117 : !Now it calculates the exponential of the eigenvectors (times the factor)
118 :
119 : !--exponential of the diagonal
120 0 : ww(:) = exp( ww(:)*factor )
121 :
122 : !--construction exponential matrix
123 0 : mat_a = zero
124 0 : mat_a(:,1) = ww(:)
125 0 : mat_a(:,:) = cshift(array=mat_a,shift=(/ (-ii,ii=0,mat_a_size) /), dim=2 )
126 :
127 :
128 : !uu.exp(ww*factor)
129 0 : mat_a(:,:) = matmul(uu,mat_a)
130 :
131 : !the inverse of the eigenvectors matrix
132 0 : call ZGETRF( mat_a_size, mat_a_size, uu,mat_a_size, ipvt, info )
133 0 : if (info/=0) then
134 0 : write(msg,'(a,i4)')'Wrong value for rwork ',info
135 0 : ABI_BUG(msg)
136 : end if
137 :
138 0 : call ZGETRI( mat_a_size, uu, mat_a_size, ipvt, work, lwork, info )
139 0 : if (info/=0) then
140 0 : write(msg,'(a,i4)')'Wrong value for rwork ',info
141 0 : ABI_BUG(msg)
142 : end if
143 :
144 : !(uu.exp(ww*factor)).uu-1
145 0 : mat_a = matmul(mat_a,uu)
146 :
147 0 : ABI_FREE(ww)
148 0 : ABI_FREE(uu)
149 :
150 0 : end subroutine exp_mat_cx
151 : !!***
152 :
153 :
154 :
155 :
156 : END MODULE m_exp_mat
157 : !!***
|
