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
zpoly.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 Matrix polynom
11 C
12 C Input/Output parameters:
13 C
14 C M - size of input square matrix (input)
15 C A - matrix (input)
16 C LDA - leading dimension of A (int)(input)
17 C N - size of coefficient array V (input)
18 C V - coefficient array (input)
19 C P - result square matrix (output) = V(1)*I + V(2)*A + V(3)*A^2 + ... + V(N)*A^(N-1)
20 C LDP - leading dimension of P (int)(input)
21 C B - working array of size NPOLY (M, N)
22 
23  SUBROUTINE zpoly (M, A, LDA, N, V, P, LDP, B)
24 CDEC$ IF DEFINED (FTN_EXPORTS)
25 CDEC$ ATTRIBUTES DLLEXPORT::ZPOLY
26 CDEC$ ENDIF
27  INTEGER m, n, lda, ldp
28  DOUBLE COMPLEX a(lda*m), v(n), p(ldp*m), b(1)
29  INTEGER i, k, ns, nr, nqsr, mm, nrm, nwrk
30  DOUBLE COMPLEX cone /(1.d0,0.d0)/, czero /(0.d0,0.d0)/
31  DOUBLE PRECISION q
32  CHARACTER trans /'N'/
33  INTEGER npoly
34  INTEGER floor, ceiling
35 
36  IF (m .LE. 0) RETURN
37  mm = m * m
38 
39  CALL zscal(mm, czero, p, 1)
40  IF (n .LE. 0) RETURN
41 
42  DO 10 i = 0, m-1
43  p(i*(ldp+1)+1) = v(1)
44 10 CONTINUE
45 
46  IF (n .EQ. 1) RETURN
47 
48  q = dfloat(n - 1) ! b(0) is V(1)
49  ns = ceiling(dsqrt(q))
50  nr = floor(q / dfloat(ns))
51  nqsr = n - 1 - ns * nr ! q - sr
52  nrm = (nr + 1) * mm + 1 ! the 1st el. of the last matrix
53 
54  nwrk = npoly(m, n)
55  DO 15 i = mm + 1, nwrk
56  b(i) = czero
57 15 CONTINUE
58 
59  CALL zcopym(m, m, a, lda, b, m) ! the first of these matrices
60  ! will contain powers of A
61  DO 20 k = 1, nr
62  DO 30 i = 1, mm, m + 1
63  b(mm * k + i) = v(ns * k + 1) ! b(i)*I
64 30 CONTINUE
65 20 CONTINUE
66 
67  DO 40 i = 1, ns - 1
68  CALL zaxpym(m, m, v(i + 1), b, m, p, ldp)
69 
70  DO 50 k = 1, nr - 1
71  CALL zaxpy(mm, v(ns * k + i + 1), b, 1,
72  1 b(k * mm + 1), 1)
73 50 CONTINUE
74 
75  IF (i .LE. nqsr) THEN
76  CALL zaxpy(mm, v(ns * nr + i + 1), b, 1,
77  1 b(nr * mm + 1), 1)
78  ENDIF
79 
80  CALL zcopy(mm, b, 1, b(nrm), 1)
81  CALL zgemm(trans, trans, m, m, m, cone, b(nrm), m, a, lda,
82  1 czero, b, m)
83 40 CONTINUE
84 
85  CALL zgemm(trans, trans, m, m, m, cone, b(mm + 1), m,
86  1 b, m, cone, p, ldp)
87 
88  IF (nr .GT. 1) CALL zcopy(mm, b, 1, b(nrm), 1)
89 
90  DO 60 k = 2, nr
91  CALL zcopy(mm, b, 1, b(mm + 1), 1)
92  CALL zgemm(trans, trans, m, m, m, cone, b(mm + 1), m,
93  1 b(nrm), m, czero, b, m)
94  CALL zgemm(trans, trans, m, m, m, cone, b(k * mm + 1), m,
95  1 b, m, cone, p, ldp)
96 60 CONTINUE
97 
98  RETURN
99  END !SUBROUTINE ZPOLY