forked from davecats/activeBarriers
-
Notifications
You must be signed in to change notification settings - Fork 2
/
lapacksvd.cpl
executable file
·85 lines (79 loc) · 2.67 KB
/
lapacksvd.cpl
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
! compile with -fp-trap-all=nodivzero with Intel
USE complex
#link "-llapack -lblas"
INLINE ARRAY(A.LO..A.HI) OF COMPLEX FUNCTION eigv(COMPLEX A(*,*))
INTEGER SDIM,info
LA=LENGTH(A)
IF LA#LENGTH(A(LO)) THEN ERROR "matrix is not square"
LWORK=3*LA
COMPLEX WORK(1..LWORK)
REAL RWORK(1..LA)
COMPLEX Acopy(0..LA-1,0..LA-1)=A(LO+*,LO+*)
FORTRANCALL zgees("N","N",NULL,LA,Acopy,STRIDEOF(Acopy),SDIM,RESULT,NULL,LA,WORK,LWORK,RWORK,NULL,info)
IF info#0 THEN ERROR "zgees: " info
END eigv
INLINE ARRAY(A.LO..A.HI) OF COMPLEX FUNCTION eigv(REAL A(*,*))
INTEGER SDIM,info
LA=LENGTH(A)
IF LA#LENGTH(A(LO)) THEN ERROR "matrix is not square"
LWORK=3*LA
REAL WORK(1..LWORK)
REAL WR(1..LA),WI(1..LA)
REAL Acopy(0..LA-1,0..LA-1)=A(LO+*,LO+*)
FORTRANCALL dgees("N","N",NULL,LA,Acopy,STRIDEOF(Acopy),SDIM,WR,WI,NULL,LA,WORK,LWORK,NULL,info)
IF info#0 THEN ERROR "dgees: " info
RESULT(*+LO-1).REAL=WR; RESULT(*+LO-1).IMAG=WI
END eigv
SUBROUTINE svd(ARRAY(*,*) OF REAL A; POINTER TO ARRAY(*) OF REAL SVD)
INTEGER info
LA1=LENGTH(A)
LA2=LENGTH(A(LO))
LWORK=LA1*LA2*50
REAL WORK(0..LWORK)
REAL Acopy(0..LA1-1,0..LA2-1)=A(LO+*,LO+*)
FORTRANCALL dgesvd("N","N",LA1,LA2,Acopy,STRIDEOF(Acopy),SVD,NULL,LA1,NULL,LA1,WORK,LWORK,info)
IF info#0 THEN ERROR "dgesvd: " info
END svd
SUBROUTINE svdv(ARRAY(*,*) OF REAL A; POINTER TO ARRAY(*) OF REAL SVD; POINTER TO ARRAY(*,*) OF REAL L,R)
INTEGER info
LA1=LENGTH(A)
LA2=LENGTH(A(LO))
LWORK=LA1*LA2*50
REAL WORK(0..LWORK)
REAL Acopy(0..LA1-1,0..LA2-1)=A(LO+*,LO+*)
FORTRANCALL dgesvd("A","A",LA1,LA2,Acopy,STRIDEOF(Acopy),SVD,L,LA1,R,LA1,WORK,LWORK,info)
IF info#0 THEN WRITE "dgesvd: " info
END svdv
SUBROUTINE eigs(ARRAY(*,*) OF REAL A; POINTER TO ARRAY(*) OF COMPLEX SVD; POINTER TO ARRAY(*,*) OF COMPLEX L,R)
INTEGER info
LA1=LENGTH(A)
LA2=LENGTH(A(LO))
LWORK=LA1*LA2*50
REAL WORK(0..LWORK)
REAL Acopy(0..LA1-1,0..LA2-1)=A(LO+*,LO+*)
REAL WR(0..LA1-1),WI(0..LA1-1)
REAL LL(0..LA1-1,0..LA1-1), RR(0..LA1-1,0..LA1-1)
FORTRANCALL dgeev("V","V",LA2,Acopy,LA1,WR,WI,LL,LA1,RR,LA1,WORK,LWORK,info)
SVD(*+LO).REAL=WR; SVD(*+LO).IMAG=WI
INTEGER i=0
LOOP WHILE i<=SVD.HI
IF ABS(SVD(i).IMAG)>1E-10 THEN
!L(i+LO,*+LO) = LL(*,i+LO) + I*LL(*,i+LO+1)
!L(i+LO+1,*+LO) = LL(*,i+LO) - I*LL(*,i+LO+1)
!R(i+LO,*+LO) = RR(*,i+LO) + I*RR(*,i+LO+1)
!R(i+LO+1,*+LO) = RR(*,i+LO) - I*RR(*,i+LO+1)
L(i+LO,*+LO) = LL(i+LO,*) - I*LL(i+LO+1,*)
L(i+LO+1,*+LO) = LL(i+LO,*) + I*LL(i+LO+1,*)
R(i+LO,*+LO) = RR(i+LO,*) - I*RR(i+LO+1,*)
R(i+LO+1,*+LO) = RR(i+LO,*) + I*RR(i+LO+1,*)
INC i
ELSE
!L(i+LO,*+LO) = LL(*,i+LO)
!R(i+LO,*+LO) = RR(*,i+LO)
L(i+LO,*+LO) = LL(i+LO,*)
R(i+LO,*+LO) = RR(i+LO,*)
END IF
INC i
REPEAT
IF info#0 THEN WRITE "dgesvd: " info
END eigs