CVM Class Library  8.1
This C++ class library encapsulates concepts of vector and different matrices including square, band, symmetric and hermitian ones in Euclidean space of real and complex numbers.
 All Classes Files Functions Variables Typedefs Friends Macros Pages
zmexp.f
Go to the documentation of this file.
1 C CVM Class Library
2 C http://cvmlib.com
3 C
4 C Copyright Sergei Nikolaev 1992-2013
5 C Distributed under the Boost Software License, Version 1.0.
6 C (See accompanying file LICENSE_1_0.txt or copy at
7 C http://www.boost.org/LICENSE_1_0.txt)
8 C
9 C
10 C Square matrix exponent
11 C
12 C Input/Output parameters:
13 C
14 C M - size of matrix (int)(input)
15 C A - matrix (complex)(input)
16 C LDA - leading dimesion of A (int)(input)
17 C EA - exponent of A (complex)(output)
18 C LDE - leading dimesion of EA (int)(input)
19 C R - working array of size NR (real)(input). IT MUST BE FILLED WITH ZEROES
20 C IR - working array of size NI (int)(input)
21 C NR - working array R size (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
22 C NI - working array IR size (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
23 C NQ - row length (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
24 C J - A measure (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
25 C ISHEM - whether matrix A is hermitian (int)(input). 1 - true, 0 - false
26 C WORK - working array of size LWORK (real)(input)
27 C LWORK - length of working array WORK (int)(input) - usually 64 * M
28 
29  SUBROUTINE zmexp (M, A, LDA, EA, LDE, R, IR, NR, NI, NQ, J,
30  1 ishem, work, lwork) ! referenced in
31  ! symmetric case only
32 CDEC$ IF DEFINED (FTN_EXPORTS)
33 CDEC$ ATTRIBUTES DLLEXPORT::ZMEXP
34 CDEC$ ENDIF
35  INTEGER m, lda, lde
36  DOUBLE COMPLEX a(lda*m), ea(lde*m)
37  INTEGER nr, ni, nq, j, lwork
38  LOGICAL ishem
39  DOUBLE COMPLEX r(nr), work(lwork)
40  INTEGER ir(ni)
41 
42  INTEGER i, iq, nqc, nb, info, nnr
43  INTEGER mm, mm1, mm12, mm13, mm14, mnb, mw
44  DOUBLE PRECISION two /2.d0/
45  DOUBLE COMPLEX czero /(0.d0, 0.d0)/, cone /(1.d0, 0.d0)/
46  DOUBLE COMPLEX tj
47  DOUBLE PRECISION c
48  CHARACTER trans /'N'/
49  CHARACTER uplo /'U'/
50  LOGICAL even
51  INTEGER floor, ceiling
52 
53  IF (m .LE. 0) RETURN
54  IF (m .EQ. 1) THEN
55  ea(1) = zexp(a(1))
56  RETURN
57  END IF
58 
59  tj = (1.d0, 0.d0)
60  c = 0.5d0
61 
62  mm = m * m
63  mm1 = mm + 1
64  mm12 = mm1 + mm
65  mm13 = mm12 + mm
66  mm14 = mm13 + mm
67 
68 C R = CZERO
69 
70  nqc = floor(dfloat(nq) / two) + 1
71  mnb = mm14 + 2 * nqc
72  nnr = floor(dfloat(nqc - 1) /
73  1 dfloat(ceiling(dsqrt(dfloat(nqc - 1)))))
74  nb = (2 * nnr + 2) * mm
75  mw = mnb + nb
76 
77  even = .true.
78  i = 1
79  r(mm14) = cone
80  r(mm14 + nqc) = dcmplx(c)
81  DO 40 iq = 2, nq
82  c = c * dfloat(nq - iq + 1) / dfloat((2 * nq - iq + 1) * iq)
83 
84  IF (even) THEN
85  r(mm14 + i) = dcmplx(c)
86  ELSE
87  r(mm14 + nqc + i) = dcmplx(c)
88  i = i + 1
89  ENDIF
90 
91  even = .NOT. even
92 40 CONTINUE
93 
94  CALL zcopym(m, m, a, lda, ea, lde)
95  IF (j .GT. 0) THEN
96  tj = dcmplx(two ** dfloat(-j))
97  CALL zscalm(m, m, tj, ea, lde)
98  ENDIF
99 
100  CALL zcopym(m, m, ea, lde, r(mm1), m)
101  CALL zgemm(trans, trans, m, m, m, cone, r(mm1), m,
102  1 ea, lde, czero, r, m)
103  CALL zscal(mm, czero, r(mm1), 1)
104 
105 C U and V matrices calculation...
106  CALL zpoly2(m, r, nqc, r(mm14), r(mm14 + nqc), r(mm1), r(mm12),
107  1 r(mnb), r(mnb + (nnr + 2) * mm))
108  CALL zcopym(m, m, r(mm1), m, ea, lde)
109  CALL zgemm(trans, trans, m, m, m, tj,
110  1 a, m, r(mm12), m, cone, ea, lde) ! N matrix
111  CALL zgemm(trans, trans, m, m, m, -tj,
112  1 a, m, r(mm12), m, cone, r(mm1), m) ! D matrix
113  CALL zcopy(mm, r(mm1), 1, r(mm13), 1) ! copy of D
114  CALL zcopym(m, m, ea, lde, r(mm12), m) ! copy of N
115 
116  IF (ishem) THEN
117  CALL zhetrf(uplo, m, r(mm1), m, ir, work, lwork, info)
118  CALL zcopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
119  CALL zhetrs(uplo, m, m, r(mm1), m, ir, ea, lde, info)
120  CALL zherfs(uplo, m, m, r(mm13), m, r, m, ir,
121  1 r(mm12), m, ea, lde, r(mw), r(mw + m),
122  2 r(mw + 2*m), work, info)
123  ELSE
124  CALL zgetrf(m, m, r(mm1), m, ir, info) ! factorize
125  CALL zcopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
126  ! F matrix (EA) is the solution of the equation DF = N
127  CALL zgetrs(trans, m, m, r(mm1), m, ir, ea, lde, info)
128  CALL zgerfs(trans, m, m, r(mm13), m, r, m, ir,
129  1 r(mm12), m, ea, lde, r(mw), r(mw + m), r(mw + 2*m),
130  2 r(mm1), info)
131  ENDIF
132 
133 
134  DO 50 i = 1, j
135  CALL zcopym(m, m, ea, lde, r, m)
136  CALL zcopym(m, m, ea, lde, r(mm1), m)
137  CALL zgemm(trans, trans, m, m, m, cone,
138  1 r, m, r(mm1), m, czero, ea, lde) ! F = F * F
139 50 CONTINUE
140 
141  RETURN
142  END !SUBROUTINE ZMEXP
143 
144 
145 C Input/Output parameters:
146 C
147 C M - size of matrix (int)(input)
148 C A - matrix (complex)(input)
149 C LDA - leading dimesion of A (int)(input)
150 C TOL - tolerance (real)(input)
151 C NR - working array R size (int)(output).
152 C NI - working array IR size (int)(output).
153 C NQ - row length (int)(output).
154 C J - A measure (int)(output).
155 
156  SUBROUTINE zmexpc (M, A, LDA, TOL, NR, NI, NQ, J)
157 CDEC$ IF DEFINED (FTN_EXPORTS)
158 CDEC$ ATTRIBUTES DLLEXPORT::ZMEXPC
159 CDEC$ ENDIF
160  INTEGER m, lda
161  DOUBLE COMPLEX a(lda*m)
162  DOUBLE PRECISION tol
163  INTEGER nr, ni, nq, j
164  DOUBLE PRECISION anrm, tmp, eq
165  DOUBLE PRECISION one /1.d0/, two /2.d0/
166  DOUBLE PRECISION twol /0.69314718055994530941723212145818d0/
167  INTEGER nb, mm, nqc
168  DOUBLE PRECISION zinfnm
169  INTEGER floor, ceiling
170  DOUBLE PRECISION tiny
171 
172  ni = m * 2
173  mm = m * m
174 
175  anrm = zinfnm(m, m, a, lda)
176  nq = 1
177  eq = 0.16666666666666666666666666666667d0 ! 1/6
178  DO 10 WHILE (.true.)
179  tmp = eq * anrm
180  IF (tmp * dexp(tmp) .LT. tol .AND. nq .GT. 1) goto 15
181  nq = nq + 1
182  eq = eq / dfloat(16 * (4 * nq * nq - 1))
183 10 CONTINUE
184 
185 15 nq = nq + 1
186 
187  j = 0
188  IF (anrm .GT. tiny(one)) THEN
189  j = 1 + ceiling(dlog(anrm) / twol)
190  ENDIF
191 
192  nqc = floor(dfloat(nq) / two) + 1
193  nb = floor(dfloat(nqc - 1) /
194  1 dfloat(ceiling(dsqrt(dfloat(nqc - 1)))))
195  nb = (2 * nb + 2) * mm
196 
197  nr = 4 * mm + 2 * nqc + nb + 5 * m
198  RETURN
199  END !SUBROUTINE ZMEXPC
200 
201 
202  SUBROUTINE zpoly2 (M, A, N, V1, V2, P1, P2, B, B2)
203  INTEGER m, n
204  DOUBLE COMPLEX a(m*m), v1(n), v2(n), p1(m*m), p2(m*m), b(1), b2(1)
205  DOUBLE COMPLEX cone /(1.d0, 0.d0)/, czero /(0.d0, 0.d0)/
206  INTEGER i, k, ns, nr, nqsr, mm, nrm
207  DOUBLE PRECISION q
208  CHARACTER trans /'N'/
209  INTEGER floor, ceiling
210 
211  mm = m * m
212  DO 10 i = 1, mm, m + 1
213  p1(i) = v1(1)
214  p2(i) = v2(1)
215 10 CONTINUE
216 
217  q = dfloat(n - 1) ! b(0) is V(1)
218  ns = ceiling(dsqrt(q))
219  nr = floor(q / dfloat(ns))
220  nqsr = n - 1 - ns * nr ! q - sr
221  nrm = (nr + 1) * mm + 1 ! the 1st el. of the last matrix
222 
223  CALL zcopy(mm, a, 1, b, 1) ! the first of these matrices
224  ! will contain powers of A
225  DO 20 k = 1, nr
226  DO 30 i = 1, mm, m + 1
227  b(mm * k + i) = v1(ns * k + 1) ! b(i)*I
228  b2(mm * (k - 1) + i) = v2(ns * k + 1) ! b(i)*I
229 30 CONTINUE
230 20 CONTINUE
231 
232  DO 40 i = 1, ns - 1
233  CALL zaxpy(mm, v1(i + 1), b, 1, p1, 1)
234  CALL zaxpy(mm, v2(i + 1), b, 1, p2, 1)
235 
236  DO 50 k = 1, nr - 1
237  CALL zaxpy(mm, v1(ns * k + i + 1), b, 1,
238  1 b(k * mm + 1), 1)
239  CALL zaxpy(mm, v2(ns * k + i + 1), b, 1,
240  1 b2((k - 1) * mm + 1), 1)
241 50 CONTINUE
242 
243  IF (i .LE. nqsr) THEN
244  CALL zaxpy(mm, v1(ns * nr + i + 1), b, 1,
245  1 b(nr * mm + 1), 1)
246  CALL zaxpy(mm, v2(ns * nr + i + 1), b, 1,
247  1 b2((nr - 1) * mm + 1), 1)
248  ENDIF
249 
250  CALL zcopy(mm, b, 1, b(nrm), 1)
251  CALL zgemm(trans, trans, m, m, m, cone, b(nrm), m, a, m,
252  1 czero, b, m)
253 40 CONTINUE
254 
255  CALL zgemm(trans, trans, m, m, m, cone, b(mm + 1), m,
256  1 b, m, cone, p1, m)
257  CALL zgemm(trans, trans, m, m, m, cone, b2, m,
258  1 b, m, cone, p2, m)
259 
260  IF (nr .GT. 1) CALL zcopy(mm, b, 1, b(nrm), 1)
261 
262  DO 60 k = 2, nr
263  CALL zcopy(mm, b, 1, b(mm + 1), 1)
264  CALL zgemm(trans, trans, m, m, m, cone, b(mm + 1), m,
265  1 b(nrm), m, czero, b, m)
266  CALL zgemm(trans, trans, m, m, m, cone,
267  1 b(k * mm + 1), m, b, m, cone, p1, m)
268  CALL zgemm(trans, trans, m, m, m, cone,
269  1 b2((k - 1) * mm + 1), m, b, m, cone, p2, m)
270 60 CONTINUE
271 
272  RETURN
273  END !SUBROUTINE ZPOLY2
Generated on Tue Jul 22 2014 18:52:07 for CVM Class Library by   doxygen 1.8.1.2