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
cmexp.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 (real)(input)
16 C LDA - leading dimesion of A (int)(input)
17 C EA - exponent of A (real)(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 cmexp (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::CMEXP
34 CDEC$ ENDIF
35  INTEGER m, lda, lde
36  COMPLEX a(lda*m), ea(lde*m)
37  INTEGER nr, ni, nq, j, lwork
38  LOGICAL ishem
39  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  REAL two /2./
45  DOUBLE PRECISION dtwo /2.d0/
46  COMPLEX czero /(0., 0.)/, cone /(1., 0.)/
47  COMPLEX tj
48  REAL c
49  CHARACTER trans /'N'/
50  CHARACTER uplo /'U'/
51  LOGICAL even
52  INTEGER floor, ceiling
53 
54  IF (m .LE. 0) RETURN
55  IF (m .EQ. 1) THEN
56  ea(1) = cexp(a(1))
57  RETURN
58  END IF
59 
60  tj = (1., 0.)
61  c = 0.5
62 
63  mm = m * m
64  mm1 = mm + 1
65  mm12 = mm1 + mm
66  mm13 = mm12 + mm
67  mm14 = mm13 + mm
68 
69  nqc = floor(dfloat(nq) / dtwo) + 1
70  mnb = mm14 + 2 * nqc
71  nnr = floor(dfloat(nqc - 1) /
72  1 dfloat(ceiling(dsqrt(dfloat(nqc - 1)))))
73  nb = (2 * nnr + 2) * mm
74  mw = mnb + nb
75 
76 
77  even = .true.
78  i = 1
79  r(mm14) = cone
80  r(mm14 + nqc) = cmplx(c)
81  DO 40 iq = 2, nq
82  c = c * float(nq - iq + 1) / float((2 * nq - iq + 1) * iq)
83 
84  IF (even) THEN
85  r(mm14 + i) = cmplx(c)
86  ELSE
87  r(mm14 + nqc + i) = cmplx(c)
88  i = i + 1
89  ENDIF
90 
91  even = .NOT. even
92 40 CONTINUE
93 
94  CALL ccopym(m, m, a, lda, ea, lde)
95  IF (j .GT. 0) THEN
96  tj = cmplx(two ** float(-j))
97  CALL cscalm(m, m, tj, ea, lde)
98  ENDIF
99 
100  CALL ccopym(m, m, ea, lde, r(mm1), m)
101  CALL cgemm(trans, trans, m, m, m, cone, r(mm1), m,
102  1 ea, lde, czero, r, m)
103  CALL cscal(mm, czero, r(mm1), 1)
104 
105 C U and V matrices calculation...
106  CALL cpoly2(m, r, nqc, r(mm14), r(mm14 + nqc), r(mm1), r(mm12),
107  1 r(mnb), r(mnb + (nnr + 2) * mm))
108  CALL ccopym(m, m, r(mm1), m, ea, lde)
109  CALL cgemm(trans, trans, m, m, m, tj,
110  1 a, lda, r(mm12), m, cone, ea, lde) ! N matrix
111  CALL cgemm(trans, trans, m, m, m, -tj,
112  1 a, lda, r(mm12), m, cone, r(mm1), m) ! D matrix
113  CALL ccopy(mm, r(mm1), 1, r(mm13), 1) ! copy of D
114  CALL ccopym(m, m, ea, lde, r(mm12), m) ! copy of N
115 
116  IF (ishem) THEN
117  CALL chetrf(uplo, m, r(mm1), m, ir, work, lwork, info)
118  CALL ccopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
119  CALL chetrs(uplo, m, m, r(mm1), m, ir, ea, lde, info)
120  CALL cherfs(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 cgetrf(m, m, r(mm1), m, ir, info) ! factorize
125  CALL ccopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
126  ! F matrix (EA) is the solution of the equation DF = N
127  CALL cgetrs(trans, m, m, r(mm1), m, ir, ea, lde, info)
128  CALL cgerfs(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  DO 50 i = 1, j
134  CALL ccopym(m, m, ea, lde, r, m)
135  CALL ccopym(m, m, ea, lde, r(mm1), m)
136  CALL cgemm(trans, trans, m, m, m, cone,
137  1 r, m, r(mm1), m, czero, ea, lde) ! F = F * F
138 50 CONTINUE
139 
140  RETURN
141  END !SUBROUTINE CMEXP
142 
143 
144 C Input/Output parameters:
145 C
146 C M - size of matrix (int)(input)
147 C A - matrix (real)(input)
148 C LDA - leading dimesion of A (int)(input)
149 C TOL - tolerance (real)(input)
150 C NR - working array R size (int)(output)
151 C NI - working array IR size (int)(output)
152 C NQ - row length (int)(output)
153 C J - A measure (int)(output)
154 
155 
156  SUBROUTINE cmexpc (M, A, LDA, TOL, NR, NI, NQ, J)
157 CDEC$ IF DEFINED (FTN_EXPORTS)
158 CDEC$ ATTRIBUTES DLLEXPORT::CMEXPC
159 CDEC$ ENDIF
160  INTEGER m, lda
161  COMPLEX a(lda*m)
162  REAL tol
163  INTEGER nr, ni, nq, j
164 
165  REAL anrm, tmp, eq
166  REAL one /1./
167  DOUBLE PRECISION dtwo /2.d0/
168  DOUBLE PRECISION dtwol /0.69314718055994530941723212145818d0/
169  INTEGER nb, mm, nqc
170  REAL cinfnm
171  INTEGER floor, ceiling
172  REAL tiny
173 
174  ni = m * 2
175  mm = m * m
176 
177  anrm = cinfnm(m, m, a, lda)
178  nq = 1
179  eq = 0.16666666666666666666666666666667 ! 1/6
180  DO 10 WHILE (.true.)
181  tmp = eq * anrm
182  IF (tmp * exp(tmp) .LT. tol .AND. nq .GT. 1) goto 15
183  nq = nq + 1
184  eq = eq / float(16 * (4 * nq * nq - 1))
185 10 CONTINUE
186 
187 15 nq = nq + 1
188 
189  j = 0
190  IF (anrm .GT. tiny(one)) THEN
191  j = 1 + ceiling(dlog(dble(anrm)) / dtwol)
192  ENDIF
193 
194  nqc = floor(dfloat(nq) / dtwo) + 1
195  nb = floor(dfloat(nqc - 1) /
196  1 dfloat(ceiling(dsqrt(dfloat(nqc - 1)))))
197  nb = (2 * nb + 2) * mm
198 
199  nr = 4 * mm + 2 * nqc + nb + 5 * m
200  RETURN
201  END !SUBROUTINE CMEXPC
202 
203 
204 
205  SUBROUTINE cpoly2 (M, A, N, V1, V2, P1, P2, B, B2)
206  INTEGER m, n
207  COMPLEX a(m*m), v1(n), v2(n), p1(m*m), p2(m*m), b(1), b2(1)
208  COMPLEX cone /(1., 0.)/, czero /(0., 0.)/
209  INTEGER i, k, ns, nr, nqsr, mm, nrm
210  REAL q
211  CHARACTER trans /'N'/
212  INTEGER floor, ceiling
213 
214  mm = m * m
215  DO 10 i = 1, mm, m + 1
216  p1(i) = v1(1)
217  p2(i) = v2(1)
218 10 CONTINUE
219 
220  q = float(n - 1) ! b(0) is V(1)
221  ns = ceiling(dsqrt(dble(q)))
222  nr = floor(dble(q) / dfloat(ns))
223  nqsr = n - 1 - ns * nr ! q - sr
224  nrm = (nr + 1) * mm + 1 ! the 1st el. of the last matrix
225 
226  CALL ccopy(mm, a, 1, b, 1) ! the first of these matrices
227  ! will contain powers of A
228  DO 20 k = 1, nr
229  DO 30 i = 1, mm, m + 1
230  b(mm * k + i) = v1(ns * k + 1) ! b(i)*I
231  b2(mm * (k - 1) + i) = v2(ns * k + 1) ! b(i)*I
232 30 CONTINUE
233 20 CONTINUE
234 
235  DO 40 i = 1, ns - 1
236  CALL caxpy(mm, v1(i + 1), b, 1, p1, 1)
237  CALL caxpy(mm, v2(i + 1), b, 1, p2, 1)
238 
239  DO 50 k = 1, nr - 1
240  CALL caxpy(mm, v1(ns * k + i + 1), b, 1,
241  1 b(k * mm + 1), 1)
242  CALL caxpy(mm, v2(ns * k + i + 1), b, 1,
243  1 b2((k - 1) * mm + 1), 1)
244 50 CONTINUE
245 
246  IF (i .LE. nqsr) THEN
247  CALL caxpy(mm, v1(ns * nr + i + 1), b, 1,
248  1 b(nr * mm + 1), 1)
249  CALL caxpy(mm, v2(ns * nr + i + 1), b, 1,
250  1 b2((nr - 1) * mm + 1), 1)
251  ENDIF
252 
253  CALL ccopy(mm, b, 1, b(nrm), 1)
254  CALL cgemm(trans, trans, m, m, m, cone, b(nrm), m, a, m,
255  1 czero, b, m)
256 40 CONTINUE
257 
258  CALL cgemm(trans, trans, m, m, m, cone, b(mm + 1), m,
259  1 b, m, cone, p1, m)
260  CALL cgemm(trans, trans, m, m, m, cone, b2, m,
261  1 b, m, cone, p2, m)
262 
263  IF (nr .GT. 1) CALL ccopy(mm, b, 1, b(nrm), 1)
264 
265  DO 60 k = 2, nr
266  CALL ccopy(mm, b, 1, b(mm + 1), 1)
267  CALL cgemm(trans, trans, m, m, m, cone, b(mm + 1), m,
268  1 b(nrm), m, czero, b, m)
269  CALL cgemm(trans, trans, m, m, m, cone,
270  1 b(k * mm + 1), m, b, m, cone, p1, m)
271  CALL cgemm(trans, trans, m, m, m, cone,
272  1 b2((k - 1) * mm + 1), m, b, m, cone, p2, m)
273 60 CONTINUE
274 
275  RETURN
276  END !SUBROUTINE CPOLY2
Generated on Tue Jul 22 2014 18:52:07 for CVM Class Library by   doxygen 1.8.1.2