TABLE OF CONTENTS
ABINIT/m_exp_mat [ Modules ]
NAME
m_exp_mat
FUNCTION
This subroutine calculate the exponential of a matrix
COPYRIGHT
Copyright (C) 2002-2024 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
SOURCE
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 !===========================================================
m_exp_mat/exp_mat_cx [ Functions ]
[ Top ] [ m_exp_mat ] [ Functions ]
NAME
exp_mat_cx
FUNCTION
Returns the exponential of a complex matrix multiplied a scalar
INPUTS
mat_a = the complex matrix mat_a = the size of mat_a factor = a real factor multiplying at the exponent
OUTPUT
exp_mat_cx = its exponential is returned in the same matrix
SOURCE
62 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(dpc),intent(inout) :: mat_a(:,:) 69 !complex(dpc) :: 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 integer :: ipvt(mat_a_size) 79 character(len=500) :: msg 80 real(dp) :: rwork(2*maxsize) 81 complex(dpc),allocatable :: ww(:),uu(:,:) 82 complex(dpc) :: work(lwork),vl(1,1) 83 ! ********************************************************************* 84 85 mat_a_size = max(1,size(mat_a,dim=1)) 86 ABI_MALLOC(ww,(mat_a_size)) 87 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 vl,1,uu, mat_a_size, work, lwork, rwork, info) 92 if (info/=0) then 93 write(msg,'(a,i4)')'Wrong value for rwork ',info 94 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 ww(:) = exp( ww(:)*factor ) 121 122 !--construction exponential matrix 123 mat_a = zero 124 mat_a(:,1) = ww(:) 125 mat_a(:,:) = cshift(array=mat_a,shift=(/ (-ii,ii=0,mat_a_size) /), dim=2 ) 126 127 128 !uu.exp(ww*factor) 129 mat_a(:,:) = matmul(uu,mat_a) 130 131 !the inverse of the eigenvectors matrix 132 call ZGETRF( mat_a_size, mat_a_size, uu,mat_a_size, ipvt, info ) 133 if (info/=0) then 134 write(msg,'(a,i4)')'Wrong value for rwork ',info 135 ABI_BUG(msg) 136 end if 137 138 call ZGETRI( mat_a_size, uu, mat_a_size, ipvt, work, lwork, info ) 139 if (info/=0) then 140 write(msg,'(a,i4)')'Wrong value for rwork ',info 141 ABI_BUG(msg) 142 end if 143 144 !(uu.exp(ww*factor)).uu-1 145 mat_a = matmul(mat_a,uu) 146 147 ABI_FREE(ww) 148 ABI_FREE(uu) 149 150 end subroutine exp_mat_cx