TABLE OF CONTENTS


ABINIT/m_exp_mat [ Modules ]

[ Top ] [ Modules ]

NAME

  m_exp_mat

FUNCTION

  This subroutine calculate the exponential of a  matrix

COPYRIGHT

 Copyright (C) 2002-2018 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.

PARENTS

CHILDREN

SOURCE

23 #if defined HAVE_CONFIG_H
24 #include "config.h"
25 #endif
26 
27 #include "abi_common.h"
28 
29 MODULE m_exp_mat
30 
31  use defs_basis
32  use m_abicore
33  use m_errors
34  use m_linalg_interfaces
35 
36  implicit none
37 
38  private  
39 
40  public :: exp_mat  ! exponential of a complex matrix   
41 
42 
43  interface exp_mat
44   module procedure exp_mat_cx
45  end interface exp_mat
46 
47 
48 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

PARENTS

CHILDREN

      zgeev,zgetrf,zgetri

SOURCE

 72  subroutine exp_mat_cx(mat_a,mat_a_size,factor)
 73 
 74 
 75 !This section has been created automatically by the script Abilint (TD).
 76 !Do not modify the following lines by hand.
 77 #undef ABI_FUNC
 78 #define ABI_FUNC 'exp_mat_cx'
 79 !End of the abilint section
 80 
 81   implicit none
 82   !Arguments ------------------------------------
 83   ! scalars
 84   real(dp),intent(in) ::  factor
 85   ! arrays
 86   complex(dpc),intent(inout) :: mat_a(:,:)
 87   !complex(dpc) :: exp_mat_cx(mat_a_size,mat_a_size)
 88 
 89   !Local ------------------------------------------
 90   ! scalars
 91   integer :: info,mat_a_size,ii
 92   integer,parameter :: maxsize=3
 93   integer,parameter :: lwork=(1+32)*maxsize 
 94 
 95   ! arrays
 96   integer :: ipvt(mat_a_size)
 97   character(len=500) :: msg
 98   real(dp) :: rwork(2*maxsize)
 99   complex(dpc),allocatable :: ww(:),uu(:,:)
100   complex(dpc) :: work(lwork),vl(1,1)
101   ! *********************************************************************
102 
103   mat_a_size = max(1,size(mat_a,dim=1))
104   ABI_ALLOCATE(ww,(mat_a_size))
105   ABI_ALLOCATE(uu,(mat_a_size,mat_a_size))
106 
107   !Now it calculates the eigenvalues and eigenvectors of the matrix
108   call ZGEEV('No left vectors','Vectors (right)',mat_a_size, mat_a, mat_a_size,ww,&
109     vl,1,uu, mat_a_size, work, lwork, rwork, info)
110   if (info/=0) then
111    write(msg,'(a,i4)')'Wrong value for rwork ',info
112    MSG_BUG(msg)
113   end if
114 
115   !!debbug
116   !    write(std_out,*)'mat_a',mat_a
117   !    write(std_out,*)'mat_a_size',mat_a_size
118 
119   !    write(std_out,*)'eigenvalues'
120   !    write(std_out,*)ww
121   !    write(std_out,*)'eigenvectors'
122   !    do ii=1,g_mat_size
123   !     write(std_out,*)'autov',ii
124   !     write(std_out,*)uu(:,ii)
125   !    end do
126   !    write(std_out,*)'optimal workl=',work(1)
127 
128   !    write(std_out,*)'info', info
129   !    write(std_out,*) '------------------'
130   !!enddebug
131 
132 
133 
134   !-----------------------------------------------------------
135   !Now it calculates the exponential of the eigenvectors (times the factor)
136 
137   !--exponential of the diagonal
138   ww(:) = exp( ww(:)*factor )
139 
140   !--construction exponential matrix
141   mat_a = zero
142   mat_a(:,1) = ww(:)
143   mat_a(:,:) = cshift(array=mat_a,shift=(/ (-ii,ii=0,mat_a_size) /), dim=2 )
144 
145 
146   !uu.exp(ww*factor)
147   mat_a(:,:) = matmul(uu,mat_a)
148 
149   !the inverse of the eigenvectors matrix
150   call ZGETRF( mat_a_size, mat_a_size, uu,mat_a_size, ipvt, info )
151   if (info/=0) then
152    write(msg,'(a,i4)')'Wrong value for rwork ',info
153    MSG_BUG(msg)
154   end if
155 
156   call ZGETRI( mat_a_size, uu, mat_a_size, ipvt, work, lwork, info )
157   if (info/=0) then
158    write(msg,'(a,i4)')'Wrong value for rwork ',info
159    MSG_BUG(msg)
160   end if
161 
162   !(uu.exp(ww*factor)).uu-1
163   mat_a = matmul(mat_a,uu)
164 
165   ABI_DEALLOCATE(ww)
166   ABI_DEALLOCATE(uu)
167 
168  end subroutine exp_mat_cx