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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
C/MEMBR ADD NAME=WPADE,SSI=0
c Copyright INRIA
subroutine wpade(ar,ai,ia,n,ear,eai,iea,alpha,w,ipvt,ierr)
c
c!purpose
c compute the pade approximants of the exponential of a complex
c matrix a. we scale a until the spectral radius of a*2**-m
c are smaler than one.
c
c!calling sequence
c
c subroutine wpade(ar,ai,ia,n,ear,eai,iea,alpha,w,ipvt,ierr)
c
c integer ia,n,iea,ipvt,ierr
c double precision ar,ai,ear,eai,alpha,w
c dimension ar(ia,n),ai(ia,n),ear(iea,n),eai(iea,n),w(*),ipvt(*)
c
c ar,ai : array containing the matrix a
c ia : the leading dimension of arrays a.
c n : the order of the matrices a,ea .
c ear,eai : the array that contains the n*n
c matrix exp(a).
c iea : the leading dimension of array ea.
c alpha : variable containing the maximun
c norm of the eigenvalues of a.
c w : workspace array of size 4*n +4*n*n
c ipvt : integer workspace of size n
c ierr : error indicator
c ierr= 0 if normal return
c =-4 if alpha is to big for any accuracy.
c
c common /dcoeff/ c, ndng
c double precision c(41)
c integer ndng
c
c c : array containing on return pade coefficients
c ndng : on first call ndng must be set to -1,on return
c contains degree of pade approximant
c
c!auxiliary routines
c wclmat coef wcerr (j. roche)
c wmmul dmcopy (blas.extension)
c wgeco wgesl (linpack.extension)
c sqrt (fortran)
c!
c
integer ia,n,iea,ipvt,ierr
double precision ar,ai,ear,eai,alpha,w
dimension ar(ia,n),ai(ia,n),ear(iea,n),eai(iea,n),w(*),ipvt(*)
c
common /dcoeff/ c, ndng
c internal variables
integer i,j,k,m,ndng,maxc,n2
double precision rcond,c,efact,two,zero,norm,one
dimension c(41)
c
data zero, one, two, maxc /0.0d+0,1.0d+0,2.0d+0,10/
n2=n*n
kr=1
ki=kr+n2
kw=ki+n2
c
if (ndng.ge.0) go to 10
c
c compute de pade's approximants type which is necessary to obtain
c machine precision
c
call coef(ierr)
if(ierr.ne.0) goto 170
10 m = 0
efact = one
if (alpha.le.1.0d+0) go to 90
do 20 i=1,maxc
m = m + 1
efact = efact*two
if (alpha.le.efact) go to 60
20 continue
ierr = -4
go to 170
30 m = m + 1
efact = efact*two
do 50 i=1,n
do 40 j=1,n
ar(i,j) = ar(i,j)/two
ai(i,j) = ai(i,j)/two
40 continue
50 continue
norm = norm/two
go to 115
c
c we find a matrix a'=a*2-m whith a spectral radius smaller than one.
c
60 do 80 i=1,n
do 70 j=1,n
ar(i,j) = ar(i,j)/efact
ai(i,j) = ai(i,j)/efact
70 continue
80 continue
90 continue
c
c
call wcerr(ar,ai,w,ia,n,ndng,m,maxc)
c
c
norm = zero
do 110 i=1,n
alpha = zero
do 100 j=1,n
alpha = alpha + abs(ar(i,j)) + abs(ai(i,j))
100 continue
if (alpha.gt.norm) norm = alpha
110 continue
c
c compute the inverse of the denominator of dpade's approximants.
c
115 continue
do 130 i=1,n
do 120 j=1,n
ear(i,j) = -ar(i,j)
eai(i,j) = -ai(i,j)
120 continue
130 continue
call wclmat(iea,n,ear,eai,w(kr),w(ki),n,w(kw),c,ndng)
c
c compute de l-u decomposition of n (-a) and the condition number
c pp
c
call wgeco(w(kr),w(ki), n, n, ipvt, rcond, w(kw),w(kw+n))
c
rcond=rcond**4
if ((rcond+one .le. one) .and. ((norm.gt.one) .and.
* (m.lt.maxc))) go to 30
c
c compute the numerator of dpade's approximants.
c
call wclmat(ia, n, ar,ai,ear,eai, iea, w(kw), c, ndng)
c
c compute the dpade's approximants by
c
c n (-a) x=n (a)
c pp pp
c
do 150 j=1,n
call wgesl(w(kr),w(ki), n, n, ipvt, ear(1,j),eai(1,j), 0)
150 continue
if (m.eq.0) go to 170
c
c remove the effects of normalization.
c
do 160 k=1,m
call wmmul(ear,eai,iea,ear,eai,iea,w(kr),w(ki),n,n,n,n)
call dmcopy(w(kr),n,ear,iea,n,n)
call dmcopy(w(ki),n,eai,iea,n,n)
160 continue
170 continue
return
end
|
