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
smexp.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 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 smexp (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::SMEXP
35 CDEC$ ENDIF
36  INTEGER m, lda, lde
37  REAL a(lda*m), ea(lde*m)
38  INTEGER nr, ni, nq, j, lwork
39  LOGICAL issymm
40  REAL 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  REAL zero /0./, one /1./, two /2./
46  DOUBLE PRECISION dtwo /2.d0/
47  REAL tj, 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) = exp(a(1))
56  RETURN
57  END IF
58 
59  tj = 1.
60  c = 0.5
61 
62  mm = m * m
63  mm1 = mm + 1
64  mm12 = mm1 + mm
65  mm13 = mm12 + mm
66  mm14 = mm13 + mm
67 
68  nqc = floor(dfloat(nq) / dtwo) + 1
69  mnb = mm14 + 2 * nqc
70  nnr = floor(dfloat(nqc - 1) /
71  1 dfloat(ceiling(dsqrt(dfloat(nqc - 1)))))
72  nb = (2 * nnr + 2) * mm
73  mw = mnb + nb
74 
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 * float(nq - iq + 1) / float((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 scopym(m, m, a, lda, ea, lde)
94  IF (j .GT. 0) THEN
95  tj = two ** float(-j)
96  CALL sscalm(m, m, tj, ea, lde)
97  ENDIF
98 
99  CALL scopym(m, m, ea, lde, r(mm1), m)
100  CALL sgemm(trans, trans, m, m, m, one, r(mm1), m,
101  1 ea, lde, zero, r, m)
102  CALL sscal(mm, zero, r(mm1), 1)
103 
104 C U and V matrices calculation...
105  CALL spoly2(m, r, nqc, r(mm14), r(mm14 + nqc), r(mm1), r(mm12),
106  1 r(mnb), r(mnb + (nnr + 2) * mm))
107  CALL scopym(m, m, r(mm1), m, ea, lde)
108  CALL sgemm(trans, trans, m, m, m, tj,
109  1 a, lda, r(mm12), m, one, ea, lde) ! N matrix
110  CALL sgemm(trans, trans, m, m, m, -tj,
111  1 a, lda, r(mm12), m, one, r(mm1), m) ! D matrix
112  CALL scopy(mm, r(mm1), 1, r(mm13), 1) ! copy of D
113  CALL scopym(m, m, ea, lde, r(mm12), m) ! copy of N
114 
115  IF (issymm) THEN
116  CALL ssytrf(uplo, m, r(mm1), m, ir, work, lwork, info)
117  CALL scopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
118  CALL ssytrs(uplo, m, m, r(mm1), m, ir, ea, lde, info)
119  CALL ssyrfs(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 sgetrf(m, m, r(mm1), m, ir, info) ! factorize
124  CALL scopy(mm, r(mm1), 1, r, 1) ! copy of LD/UD
125  ! F matrix (EA) is the solution of the equation DF = N
126  CALL sgetrs(trans, m, m, r(mm1), m, ir, ea, lde, info)
127  CALL sgerfs(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 scopym(m, m, ea, lde, r, m)
134  CALL scopym(m, m, ea, lde, r(mm1), m)
135  CALL sgemm(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 SMEXP
141 
142 
143 C Input/Output parameters:
144 C
145 C M - size of matrix (int)(input)
146 C A - matrix (real)(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  SUBROUTINE smexpc (M, A, LDA, TOL, NR, NI, NQ, J)
155 CDEC$ IF DEFINED (FTN_EXPORTS)
156 CDEC$ ATTRIBUTES DLLEXPORT::SMEXPC
157 CDEC$ ENDIF
158  INTEGER m, lda
159  REAL a(lda*m), tol
160  INTEGER nr, ni, nq, j
161 
162  REAL anrm, tmp, eq
163  REAL one /1./
164  DOUBLE PRECISION dtwo /2.d0/
165  DOUBLE PRECISION dtwol /0.69314718055994530941723212145818d0/
166  INTEGER nb, mm, nqc
167  REAL sinfnm
168  INTEGER floor, ceiling
169  REAL tiny
170 
171  ni = m * 2
172  mm = m * m
173 
174  anrm = sinfnm(m, m, a, lda)
175  nq = 1
176  eq = 0.16666666666666666666666666666667 ! 1/6
177  DO 10 WHILE (.true.)
178  tmp = eq * anrm
179  IF (tmp * exp(tmp) .LT. tol .AND. nq .GT. 1) goto 15
180  nq = nq + 1
181  eq = eq / float(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(dble(anrm)) / dtwol) ! G77 bug fix
189  ENDIF
190 
191  nqc = floor(dfloat(nq) / dtwo) + 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 SMEXPC
199 
200 
201  SUBROUTINE spoly2 (M, A, N, V1, V2, P1, P2, B, B2)
202  INTEGER m, n
203  REAL a(m*m), v1(n), v2(n), p1(m*m), p2(m*m)
204  REAL b(1), b2(1)
205  INTEGER i, k, ns, nr, nqsr, mm, nrm
206  REAL q, one /1./, zero /0./
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 = float(n - 1) ! b(0) is V(1)
217  ns = ceiling(dsqrt(dble(q)))
218  nr = floor(dble(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 scopy(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 saxpy(mm, v1(i + 1), b, 1, p1, 1)
233  CALL saxpy(mm, v2(i + 1), b, 1, p2, 1)
234 
235  DO 50 k = 1, nr - 1
236  CALL saxpy(mm, v1(ns * k + i + 1), b, 1,
237  1 b(k * mm + 1), 1)
238  CALL saxpy(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 saxpy(mm, v1(ns * nr + i + 1), b, 1,
244  1 b(nr * mm + 1), 1)
245  CALL saxpy(mm, v2(ns * nr + i + 1), b, 1,
246  1 b2((nr - 1) * mm + 1), 1)
247  ENDIF
248 
249  CALL scopy(mm, b, 1, b(nrm), 1)
250  CALL sgemm(trans, trans, m, m, m, one, b(nrm), m, a, m,
251  1 zero, b, m)
252 40 CONTINUE
253 
254  CALL sgemm(trans, trans, m, m, m, one, b(mm + 1), m,
255  1 b, m, one, p1, m)
256  CALL sgemm(trans, trans, m, m, m, one, b2, m,
257  1 b, m, one, p2, m)
258 
259  IF (nr .GT. 1) CALL scopy(mm, b, 1, b(nrm), 1)
260 
261  DO 60 k = 2, nr
262  CALL scopy(mm, b, 1, b(mm + 1), 1)
263  CALL sgemm(trans, trans, m, m, m, one, b(mm + 1), m,
264  1 b(nrm), m, zero, b, m)
265  CALL sgemm(trans, trans, m, m, m, one,
266  1 b(k * mm + 1), m, b, m, one, p1, m)
267  CALL sgemm(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 SPOLY2
Generated on Tue Jul 22 2014 18:52:07 for CVM Class Library by   doxygen 1.8.1.2