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
dmexp.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 TOL - tolerance (real)(input)
20 C R - working array of size NR (real)(input). IT MUST BE FILLED WITH ZEROES
21 C IR - working array of size NI (int)(input)
22 C NR - working array R size (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
23 C NI - working array IR size (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
24 C NQ - row length (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
25 C J - A measure (int)(input). IT MUST BE CALCULATED BY 'DMEXPC'
26 C ISSYMM - whether matrix A is symmetric (int)(input). 1 - true, 0 - false
27 C WORK - working array of size LWORK (real)(input)
28 C LWORK - length of working array WORK (int)(input) - usually 64 * M
29 
30  SUBROUTINE dmexp (M, A, LDA, EA, LDE, R, IR, NR, NI, NQ, J,
31  1 issymm, work, lwork) ! referenced in
32  ! symmetric case only
33 CDEC$ IF DEFINED (FTN_EXPORTS)
34 CDEC$ ATTRIBUTES DLLEXPORT::DMEXP
35 CDEC$ ENDIF
36  INTEGER m, lda, lde
37  DOUBLE PRECISION a(lda*m), ea(lde*m)
38  INTEGER nr, ni, nq, j, lwork
39  LOGICAL issymm
40  DOUBLE PRECISION r(nr), work(lwork)
41  INTEGER ir(ni)
42 
43  INTEGER i, iq, nqc, nb, info, nnr
44  INTEGER mm, mm1, mm12, mm13, mm14, mnb, mw
45  DOUBLE PRECISION zero /0.d0/, one /1.d0/, two /2.d0/
46  DOUBLE PRECISION tj, c
47  CHARACTER trans /'N'/
48  CHARACTER uplo /'U'/
49  LOGICAL even
50  INTEGER floor, ceiling
51 
52  IF (m .LE. 0) RETURN
53  IF (m .EQ. 1) THEN
54  ea(1) = dexp(a(1))
55  RETURN
56  END IF
57 
58  tj = 1.d0
59  c = 0.5d0
60 
61  mm = m * m
62  mm1 = mm + 1
63  mm12 = mm1 + mm
64  mm13 = mm12 + mm
65  mm14 = mm13 + mm
66 
67 C R = ZERO
68 
69  nqc = floor(dfloat(nq) / two) + 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  even = .true.
77  i = 1
78  r(mm14) = one
79  r(mm14 + nqc) = c
80  DO 40 iq = 2, nq
81  c = c * dfloat(nq - iq + 1) / dfloat((2 * nq - iq + 1) * iq)
82 
83  IF (even) THEN
84  r(mm14 + i) = c
85  ELSE
86  r(mm14 + nqc + i) = c
87  i = i + 1
88  ENDIF
89 
90  even = .NOT. even
91 40 CONTINUE
92 
93  CALL dcopym(m, m, a, lda, ea, lde)
94  IF (j .GT. 0) THEN
95  tj = two ** dfloat(-j)
96  CALL dscalm(m, m, tj, ea, lde)
97  ENDIF
98 
99  CALL dcopym(m, m, ea, lde, r(mm1), m)
100  CALL dgemm(trans, trans, m, m, m, one, r(mm1), m,
101  1 ea, lde, zero, r, m)
102  CALL dscal(mm, zero, r(mm1), 1)
103 
104 C U and V matrices calculation...
105  CALL dpoly2(m, r, nqc, r(mm14), r(mm14 + nqc), r(mm1), r(mm12),
106  1 r(mnb), r(mnb + (nnr + 2) * mm))
107  CALL dcopym(m, m, r(mm1), m, ea, lde)
108  CALL dgemm(trans, trans, m, m, m, tj,
109  1 a, lda, r(mm12), m, one, ea, lde) ! N matrix
110  CALL dgemm(trans, trans, m, m, m, -tj,
111  1 a, lda, r(mm12), m, one, r(mm1), m) ! D matrix
112  CALL dcopy(mm, r(mm1), 1, r(mm13), 1) ! copy of D
113  CALL dcopym(m, m, ea, lde, r(mm12), m) ! copy of N
114 
115  IF (issymm) THEN
116  CALL dsytrf(uplo, m, r(mm1), m, ir, work, lwork, info)
117  CALL dcopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
118  CALL dsytrs(uplo, m, m, r(mm1), m, ir, ea, lde, info)
119  CALL dsyrfs(uplo, m, m, r(mm13), m, r, m, ir,
120  1 r(mm12), m, ea, lde, r(mw), r(mw + m),
121  2 r(mw + 2*m), ir(m + 1), info)
122  ELSE
123  CALL dgetrf(m, m, r(mm1), m, ir, info) ! factorize
124  CALL dcopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
125  ! F matrix (EA) is the solution of the equation DF = N
126  CALL dgetrs(trans, m, m, r(mm1), m, ir, ea, lde, info)
127  CALL dgerfs(trans, m, m, r(mm13), m, r, m, ir,
128  1 r(mm12), m, ea, lde, r(mw), r(mw + m),
129  2 r(mw + 2*m), ir(m + 1), info)
130  ENDIF
131 
132  DO 50 i = 1, j
133  CALL dcopym(m, m, ea, lde, r, m)
134  CALL dcopym(m, m, ea, lde, r(mm1), m)
135  CALL dgemm(trans, trans, m, m, m, one,
136  1 r, m, r(mm1), m, zero, ea, lde) ! F = F * F
137 50 CONTINUE
138 
139  RETURN
140  END !SUBROUTINE DMEXP
141 
142 
143 C Input/Output parameters:
144 C
145 C M - size of matrix (int)(input)
146 C A - matrix (complex)(input)
147 C LDA - leading dimesion of A (int)(input)
148 C TOL - tolerance (real)(input)
149 C NR - working array R size (int)(output)
150 C NI - working array IR size (int)(output)
151 C NQ - row length (int)(output)
152 C J - A measure (int)(output)
153 
154 
155  SUBROUTINE dmexpc (M, A, LDA, TOL, NR, NI, NQ, J)
156 CDEC$ IF DEFINED (FTN_EXPORTS)
157 CDEC$ ATTRIBUTES DLLEXPORT::DMEXPC
158 CDEC$ ENDIF
159  INTEGER m, lda
160  DOUBLE PRECISION a(lda*m), tol
161  INTEGER nr, ni, nq, j
162 
163  DOUBLE PRECISION anrm, tmp, eq
164  DOUBLE PRECISION one /1.d0/, two /2.d0/
165  DOUBLE PRECISION twol /0.69314718055994530941723212145818d0/
166  INTEGER nb, mm, nqc
167  DOUBLE PRECISION dinfnm
168  INTEGER floor, ceiling
169  DOUBLE PRECISION tiny
170 
171  ni = m * 2
172  mm = m * m
173 
174  anrm = dinfnm(m, m, a, lda)
175  nq = 1
176  eq = 0.16666666666666666666666666666667d0 ! 1/6
177  DO 10 WHILE (.true.)
178  tmp = eq * anrm
179  IF (tmp * dexp(tmp) .LT. tol .AND. nq .GT. 1) goto 15
180  nq = nq + 1
181  eq = eq / dfloat(16 * (4 * nq * nq - 1))
182 10 CONTINUE
183 
184 15 nq = nq + 1
185 
186  j = 0
187  IF (anrm .GT. tiny(one)) THEN
188  j = 1 + ceiling(dlog(anrm) / twol)
189  ENDIF
190 
191  nqc = floor(dfloat(nq) / two) + 1
192  nb = floor(dfloat(nqc - 1) /
193  1 dfloat(ceiling(dsqrt(dfloat(nqc - 1)))))
194  nb = (2 * nb + 2) * mm
195 
196  nr = 4 * mm + 2 * nqc + nb + 5 * m
197  RETURN
198  END !SUBROUTINE DMEXPC
199 
200 
201  SUBROUTINE dpoly2 (M, A, N, V1, V2, P1, P2, B, B2)
202  INTEGER m, n
203  DOUBLE PRECISION a(m*m), v1(n), v2(n), p1(m*m), p2(m*m)
204  DOUBLE PRECISION b(1), b2(1)
205  INTEGER i, k, ns, nr, nqsr, mm, nrm
206  DOUBLE PRECISION q, one /1.d0/, zero /0.d0/
207  CHARACTER trans /'N'/
208  INTEGER floor, ceiling
209 
210  mm = m * m
211  DO 10 i = 1, mm, m + 1
212  p1(i) = v1(1)
213  p2(i) = v2(1)
214 10 CONTINUE
215 
216  q = dfloat(n - 1) ! b(0) is V(1)
217  ns = ceiling(dsqrt(q))
218  nr = floor(q / dfloat(ns))
219  nqsr = n - 1 - ns * nr ! q - sr
220  nrm = (nr + 1) * mm + 1 ! the 1st el. of the last matrix
221 
222  CALL dcopy(mm, a, 1, b, 1) ! the first of these matrices
223  ! will contain powers of A
224  DO 20 k = 1, nr
225  DO 30 i = 1, mm, m + 1
226  b(mm * k + i) = v1(ns * k + 1) ! b(i)*I
227  b2(mm * (k - 1) + i) = v2(ns * k + 1) ! b(i)*I
228 30 CONTINUE
229 20 CONTINUE
230 
231  DO 40 i = 1, ns - 1
232  CALL daxpy(mm, v1(i + 1), b, 1, p1, 1)
233  CALL daxpy(mm, v2(i + 1), b, 1, p2, 1)
234 
235  DO 50 k = 1, nr - 1
236  CALL daxpy(mm, v1(ns * k + i + 1), b, 1,
237  1 b(k * mm + 1), 1)
238  CALL daxpy(mm, v2(ns * k + i + 1), b, 1,
239  1 b2((k - 1) * mm + 1), 1)
240 50 CONTINUE
241 
242  IF (i .LE. nqsr) THEN
243  CALL daxpy(mm, v1(ns * nr + i + 1), b, 1,
244  1 b(nr * mm + 1), 1)
245  CALL daxpy(mm, v2(ns * nr + i + 1), b, 1,
246  1 b2((nr - 1) * mm + 1), 1)
247  ENDIF
248 
249  CALL dcopy(mm, b, 1, b(nrm), 1)
250  CALL dgemm(trans, trans, m, m, m, one, b(nrm), m, a, m,
251  1 zero, b, m)
252 40 CONTINUE
253 
254  CALL dgemm(trans, trans, m, m, m, one, b(mm + 1), m,
255  1 b, m, one, p1, m)
256  CALL dgemm(trans, trans, m, m, m, one, b2, m,
257  1 b, m, one, p2, m)
258 
259  IF (nr .GT. 1) CALL dcopy(mm, b, 1, b(nrm), 1)
260 
261  DO 60 k = 2, nr
262  CALL dcopy(mm, b, 1, b(mm + 1), 1)
263  CALL dgemm(trans, trans, m, m, m, one, b(mm + 1), m,
264  1 b(nrm), m, zero, b, m)
265  CALL dgemm(trans, trans, m, m, m, one,
266  1 b(k * mm + 1), m, b, m, one, p1, m)
267  CALL dgemm(trans, trans, m, m, m, one,
268  1 b2((k - 1) * mm + 1), m, b, m, one, p2, m)
269 60 CONTINUE
270 
271  RETURN
272  END !SUBROUTINE DPOLY2
Generated on Tue Jul 22 2014 18:52:07 for CVM Class Library by   doxygen 1.8.1.2