Branch data 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 : : !!***
|
