dpoly.f

Go to the documentation of this file.

C                  CVM Class Library
C                  http://cvmlib.com
C
C          Copyright Sergei Nikolaev 1992-2013
C Distributed under the Boost Software License, Version 1.0.
C    (See accompanying file LICENSE_1_0.txt or copy at
C          http://www.boost.org/LICENSE_1_0.txt)
C
C
C     Matrix polynom
C
C     Input/Output parameters:
C
C     M     - size of input square matrix (input)
C     A     - matrix (input)
C     LDA   - leading dimension of A (int)(input)
C     N     - size of coefficient array V (input)
C     V     - coefficient array (input)
C     P     - result square matrix (output) = V(1)*I + V(2)*A + V(3)*A^2 + ... + V(N)*A^(N-1)
C     LDP   - leading dimension of P (int)(input)
C     B     - working array of size NPOLY (M, N)

      SUBROUTINE dpoly (M, A, LDA, N, V, P, LDP, B)
CDEC$ IF DEFINED (FTN_EXPORTS)
CDEC$     ATTRIBUTES DLLEXPORT::DPOLY
CDEC$ ENDIF
      INTEGER m, n, lda, ldp
      DOUBLE PRECISION a(lda*m), v(n), p(ldp*m), b(1)
      INTEGER i, k, ns, nr, nqsr, mm, nrm, nwrk
      DOUBLE PRECISION q, one /1.d0/, zero /0.d0/
      CHARACTER trans /'N'/
      INTEGER npoly
      INTEGER floor, ceiling

      IF (m .LE. 0) RETURN
      mm = m * m

      CALL dscal(mm, zero, p, 1)
      IF (n .LE. 0) RETURN

      DO 10 i = 0, m-1
          p(i*(ldp+1)+1) = v(1)
10    CONTINUE

      IF (n .EQ. 1) RETURN

      q    = dfloat(n - 1)        ! b(0) is V(1)
      ns   = ceiling(dsqrt(q))
      nr   = floor(q / dfloat(ns))
      nqsr = n - 1 - ns * nr        ! q - sr
      nrm  = (nr + 1) * mm + 1      ! the 1st el. of the last matrix

      nwrk = npoly(m, n)
      DO 15 i = mm + 1, nwrk
          b(i) = zero
15    CONTINUE

      CALL dcopym(m, m, a, lda, b, m)   ! the first of these matrices 
                                         ! will contain powers of A
      DO 20 k = 1, nr
          DO 30 i = 1, mm, m + 1
              b(mm * k + i) = v(ns * k + 1)  ! b(i)*I
30        CONTINUE
20    CONTINUE

      DO 40 i = 1, ns - 1
          CALL daxpym(m, m, v(i + 1), b, m, p, ldp)

          DO 50 k = 1, nr - 1
              CALL daxpy(mm, v(ns * k + i + 1), b, 1, 
     1                                           b(k * mm + 1), 1)
50        CONTINUE

          IF (i .LE. nqsr) THEN
              CALL daxpy(mm, v(ns * nr + i + 1), b, 1, 
     1                                           b(nr * mm + 1), 1)
          ENDIF

          CALL dcopy(mm, b, 1, b(nrm), 1)
          CALL dgemm(trans, trans, m, m, m, one, b(nrm), m, a, lda,
     1                                           zero, b, m)
40    CONTINUE

      CALL dgemm(trans, trans, m, m, m, one, b(mm + 1), m,
     1                                           b, m, one, p, ldp)


      IF (nr .GT. 1) CALL dcopy(mm, b, 1, b(nrm), 1)

      DO 60 k = 2, nr
          CALL dcopy(mm, b, 1, b(mm + 1), 1)
          CALL dgemm(trans, trans, m, m, m, one, b(mm + 1), m, 
     1                                           b(nrm), m, zero, b, m)
          CALL dgemm(trans, trans, m, m, m, one, b(k * mm + 1), m,
     1                                           b, m, one, p, ldp)
60    CONTINUE

      RETURN
      END !SUBROUTINE DPOLY