Fortran 常用函数及操作
因为自己在学习中最常使用的是Fortran，有时候要用到mkl库中的一些函数，但是这些函数的调用参数有很多，所以将自己常用的一些记录下来，以后忘了可以
快速的查阅。
1. cheevd复厄密矩阵对角化
厄密矩阵本征值与本征矢的求解
2.getrf,getri 矩阵求逆
getrf 对一个矩阵进行LU分解
getri 计算由LU分解后矩阵的逆
1 | program main |
结果如下
1 | a = |
3.?geev一般矩阵对角化
1 | call sgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info) |
矩阵直积
最近遇到一些要批量计算的问题,其中涉及到许多Pauli矩阵的直积,这里就整理一下如何利用Fortran来计算矩阵直积1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116! Author:YuXuanLi
! E-Mail:yxli406@gmail.com
module pub
implicit none
integer xn,yn,len2,N
parameter(xn = 30,yn = 30,N = xn*yn*4,len2 = xn*yn)
complex,parameter::im = (0.0,1.0)
real,parameter::pi = 3.14159265359
complex Ham(N,N)
integer bry(8,len2)
complex sx(2,2),sy(2,2),sz(2,2),s0(2,2)k
integer info
end module pub
!========== PROGRAM START ==========================
program sol
use pub
call test()
stop
end program sol
!===================================================
Module Kron_Mod
Implicit None
contains
Subroutine Kron( A , B , H )
complex, Intent( IN ) :: A(:,:) , B(:,:)
complex, Intent( OUT ) :: H(:,:)
Integer :: i , j , m , n , p , q
m = size( A , dim = 1 )
n = size( A , dim = 2 )
p = size( B , dim = 1 )
q = size( B , dim = 2 )
Do i = 1 , m
Do j = 1 , n
H( p*(i-1)+1 : p*i , q*(j-1)+1 : q*j ) = B * A(i,j)
End Do
End Do
End Subroutine Kron
End Module Kron_Mod
!=========================================================
! Build module in here, with Block Matrix operations.
Module Kronecker
contains
! Takes in Matrices A(i,j),B(k,l), assumed 2D, returns Kronecker Product C(i*k,j*l)
function KronProd(A,B) result(C)
IMPLICIT NONE
real, dimension (:,:), intent(in) :: A, B
real, dimension (:,:), allocatable :: C
integer :: i = 0, j = 0, k = 0, l = 0
integer :: m = 0, n = 0, p = 0, q = 0
allocate(C(size(A,1)*size(B,1),size(A,2)*size(B,2)))
C = 0
do i = 1,size(A,1)
do j = 1,size(A,2)
n=(i-1)*size(B,1) + 1
m=n+size(B,1) - 1
p=(j-1)*size(B,2) + 1
q=p+size(B,2) - 1
C(n:m,p:q) = A(i,j)*B
end do
end do
end function KronProd
!----------------------------------------------------------------
! Takes in Matrices A(i,j),B(k,l), assumed 2D, returns Direct sum
! C(i+k,j+l)
function DirSum(A,B) result(C)
real, dimension (:,:), intent(in) :: A, B
real, dimension (:,:), allocatable :: C
integer :: p = 0, q = 0
allocate(C(size(A,1)+size(B,1),size(A,2)+size(B,2)))
C = 0
p = size(A,1) + size(B,1)
q = size(A,2) + size(B,2)
C(1:size(A,1),1:size(A,2)) = A
C(size(A,1)+1:p,size(A,2)+1:q) = B
return
end function DirSum
!--------------------------------------------------------
! Takes 2 vectors, A(i),B(j), returns Direct Sum C(i+j)
function VecDirSum(A,B) result(C)
real, dimension (:), intent(in) :: A, B
real, dimension (:), allocatable :: C
allocate(C(size(A)+size(B)))
C = 0
C(1:size(A)) = A
C(size(A)+1:size(A)+size(B)) = B
return
end function VecDirSUm
end module Kronecker
!===================================================
subroutine test()
use pub
use Kron_Mod
complex H1(2**2,2**2),H2(2**2,2**2)
call Pauli()
call Kron(sx,sx,H)
H2 = KronProd(sx,sx)
Do i = 1 , 2**2
Write(*,'(4(f5.1,1x))') H1( :, i)
end do
end subroutine test
!===========================================================
subroutine Pauli()
use pub
sx(1,2) = 1
sx(2,1) = 1
!----
sy(1,2) = -im
sy(2,1) = im
!-----
sz(1,1) = 1
sz(2,2) = -1
!-----
s0(1,1) = 1
s0(2,2) = 1
end subroutine Pauli
这里提供了两个模块,里面都有计算矩阵直积的函数,但有点小问题,那就是只能进行实数矩阵的直积,复数矩阵的直积暂时还没有解决.
{:.warning}
这里时网上的另外一种计算程序1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40Module Kron_Mod
Implicit None
Integer , parameter , private :: DP = Selected_Real_Kind( 9 )
contains
Subroutine Kron( A , B , H )
Real( Kind = DP ) , Intent( IN ) :: A(:,:) , B(:,:)
Real( Kind = DP ) , Intent( OUT ) :: H(:,:)
Integer :: i , j , m , n , p , q
m = size( A , dim = 1 )
n = size( A , dim = 2 )
p = size( B , dim = 1 )
q = size( B , dim = 2 )
Do i = 1 , m
Do j = 1 , n
H( p*(i-1)+1 : p*i , q*(j-1)+1 : q*j ) = B * A(i,j)
End Do
End Do
End Subroutine Kron
End Module Kron_Mod
!========================================================
Program www_fcode_cn
use Kron_Mod
Implicit None
Integer , parameter :: DP = Selected_Real_Kind( 9 )
Integer , parameter :: m=2 , n=3 , p=4, q=5 , index1 = m*p , index2 = n*q
Real(kind=DP) :: A(m,n) , B(p,q) , H(index1,index2),sx(m,m),H2(m**2,m**2)
integer :: i
sx(1,2) = 1
sx(2,1) = 1
A = reshape( (/1,2,3,4,5,6/) , (/2,3/) )
B = reshape( (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20/) , (/4,5/) )
call Kron( A , B , H )
call Kron( sx , sx , H2 )
Do i = 1 , m**2
! Write(*,'(8(f5.1,1x))') H( :, i)
Write(*,'(4(f5.1,1x))') H2( :, i)
End do
End Program www_fcode_cn
矩阵求逆
1 | subroutine inv(ndim,Amat) |
矩阵相乘
1 | ! performs matrix-matrix multiply |
读取文件行数
通常在读取文件的时候, 并不会指导文件到底有多少行, 这个子过程就是用来确定在读取到文件结尾的时候终止循环, 从而指导一共有多少行数据1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16subroutine main2()
! 读取不明行数的文件
implicit none
integer count,stat
real h1,h2,h3
open(1,file = "wavenorm.dat")
do while (.true.)
count = count + 1
read(1,*,iostat = STAT)h1,h2,h3
if(stat .ne. 0) exit ! 当这个参数不为零的时候,证明读取到文件结尾
end do
write(*,*)h1,h2,h3
write(*,*)count
close(1)
return
end subroutine
读取不明行数的文件
1 | program main |
整数/浮点数转换为字符串
1 | subroutine plot(m3) |
公众号
相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
{:.info}
|
yxliphy@gmail.com |
wechat
alipay
