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
dpoly.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 dpoly (M, A, LDA, N, V, P, LDP, B)
24 CDEC$ IF DEFINED (FTN_EXPORTS)
25 CDEC$ ATTRIBUTES DLLEXPORT::DPOLY
26 CDEC$ ENDIF
27  INTEGER m, n, lda, ldp
28  DOUBLE PRECISION a(lda*m), v(n), p(ldp*m), b(1)
29  INTEGER i, k, ns, nr, nqsr, mm, nrm, nwrk
30  DOUBLE PRECISION q, one /1.d0/, zero /0.d0/
31  CHARACTER trans /'N'/
32  INTEGER npoly
33  INTEGER floor, ceiling
34 
35  IF (m .LE. 0) RETURN
36  mm = m * m
37 
38  CALL dscal(mm, zero, p, 1)
39  IF (n .LE. 0) RETURN
40 
41  DO 10 i = 0, m-1
42  p(i*(ldp+1)+1) = v(1)
43 10 CONTINUE
44 
45  IF (n .EQ. 1) RETURN
46 
47  q = dfloat(n - 1) ! b(0) is V(1)
48  ns = ceiling(dsqrt(q))
49  nr = floor(q / dfloat(ns))
50  nqsr = n - 1 - ns * nr ! q - sr
51  nrm = (nr + 1) * mm + 1 ! the 1st el. of the last matrix
52 
53  nwrk = npoly(m, n)
54  DO 15 i = mm + 1, nwrk
55  b(i) = zero
56 15 CONTINUE
57 
58  CALL dcopym(m, m, a, lda, b, m) ! the first of these matrices
59  ! will contain powers of A
60  DO 20 k = 1, nr
61  DO 30 i = 1, mm, m + 1
62  b(mm * k + i) = v(ns * k + 1) ! b(i)*I
63 30 CONTINUE
64 20 CONTINUE
65 
66  DO 40 i = 1, ns - 1
67  CALL daxpym(m, m, v(i + 1), b, m, p, ldp)
68 
69  DO 50 k = 1, nr - 1
70  CALL daxpy(mm, v(ns * k + i + 1), b, 1,
71  1 b(k * mm + 1), 1)
72 50 CONTINUE
73 
74  IF (i .LE. nqsr) THEN
75  CALL daxpy(mm, v(ns * nr + i + 1), b, 1,
76  1 b(nr * mm + 1), 1)
77  ENDIF
78 
79  CALL dcopy(mm, b, 1, b(nrm), 1)
80  CALL dgemm(trans, trans, m, m, m, one, b(nrm), m, a, lda,
81  1 zero, b, m)
82 40 CONTINUE
83 
84  CALL dgemm(trans, trans, m, m, m, one, b(mm + 1), m,
85  1 b, m, one, p, ldp)
86 
87 
88  IF (nr .GT. 1) CALL dcopy(mm, b, 1, b(nrm), 1)
89 
90  DO 60 k = 2, nr
91  CALL dcopy(mm, b, 1, b(mm + 1), 1)
92  CALL dgemm(trans, trans, m, m, m, one, b(mm + 1), m,
93  1 b(nrm), m, zero, b, m)
94  CALL dgemm(trans, trans, m, m, m, one, b(k * mm + 1), m,
95  1 b, m, one, p, ldp)
96 60 CONTINUE
97 
98  RETURN
99  END !SUBROUTINE DPOLY