matrix.f90

CLM集合卡曼滤波数据同化算法

      subroutine TRISPE(n, d, e, x, j)
      implicit none
      integer n, j
      real d(n), e(n-1), x(n)
      integer i

      if ( j .GT. 1) x(j) = x(j)-e(j-1)*x(j-1)
      if ( j .LT. n) x(j+1) = x(j+1) - e(j)*x(j)
      do i = j+2, n
         x(i) = -e(i-1)*x(i-1)
      enddo
      
      x(n) = x(n)/d(n)
      do i=n-1, j, -1
         x(i) = x(i)/d(i) - e(i)*x(i+1)
      enddo
      
      return
      end
!******************************************************************
!*     This program requires 2n^3 flops.

!*     SSYTRD reduces a real symmetric matrix A to real symmetric
!*     tridiagonal form T by an orthogonal similarity transformation:
!*     Q**T * A * Q = T, where Q is is represented as a product of 
!*     elementary reflectors.
!*
!*     It takes: 4/3 n^3 operations.
!*
!*     If UPLO = 'L', the matrix Q is represented as a product of elementary
!*     reflectors
!*
!*     Q = H(1) H(2) . . . H(n-1).
!*
!*     Each H(i) has the form
!*
!*     H(i) = I - tau * v * v'
!*
!*     where tau is a real scalar, and v is a real vector with
!*     v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
!*     and tau in TAU(i).
!*
!*     The contents of A on exit are illustrated by the following examples
!*     with n = 5:
!*
!*     if UPLO = 'U':                       if UPLO = 'L':
!*
!*     (  d   e   v2  v3  v4 )              (  d                  )
!*     (      d   e   v3  v4 )              (  e   d              )
!*     (          d   e   v4 )              (  v1  e   d          )
!*     (              d   e  )              (  v1  v2  e   d      )
!*     (                  d  )              (  v1  v2  v3  e   d  )
!*
!*     where d and e denote diagonal and off-diagonal elements of T, and vi
!*     denotes an element of the vector defining H(i).

      subroutine APPLIQ(n,A,tau,S,v,w)
      implicit none
      integer n
      real A(n,n), tau(n-1), S(n,n), v(n), w(n)
      integer i, j, k
      real vtw

      do k= n-2, 1, -1
!         print*, k

!     Redefine the vector v, as sqrt(tau)*v, H = I - v*v'

         v(k+1) = sqrt(tau(k))
         do i=k+2, n
            v(i) = A(i,k)*v(k+1)
         enddo

!     Calculate the lower left corner, dimension is: (n-k) x k
!     Total operations:  4 (n-k) k

         do j=1, k
            w(j) = 0.0
            do i=k+1, n
               w(j) = w(j) + v(i)*S(i,j)
            enddo
         enddo

         do j=1, k
            do i=k+1, n
               S(i,j) = S(i,j) - v(i)*w(j)
            enddo
         enddo

!     Calculate the lower right trailing block, only the lower
!     triangular part, domension (n-k) x (n-k)
!     Total operations: 4 (n-k)^2

!     Operations:  2 (n-k)^2
         do i=k+1, n
            w(i) = 0.0
            do j=k+1, i
               w(i) = w(i) + S(i,j)*v(j)
            enddo
            do j=i+1, n
               w(i) = w(i) + S(j,i)*v(j)
            enddo
         enddo

         vtw = 0.0
         do i=k+1, n
            vtw = vtw + v(i)*w(i)
         enddo
         
         vtw = 0.5*vtw
         do i=k+1,n
            w(i) = w(i) - vtw*v(i)
         enddo

!     Operations: 4/2 (n-k)^2 = 2(n-k)^2

         do j= k+1, n
            do i=j, n
               S(i,j) = S(i,j) - w(i)*v(j)-v(i)*w(j)
            enddo
         enddo
               
      enddo
      return
      end
!**************************************************

      INTEGER          FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,&
                      N4 )

      CHARACTER*( * )    NAME, OPTS
      INTEGER            ISPEC, N1, N2, N3, N4

      LOGICAL            CNAME, SNAME
      CHARACTER*1        C1
      CHARACTER*2        C2, C4
      CHARACTER*3        C3
      CHARACTER*6        SUBNAM
      INTEGER            I, IC, IZ, NB, NBMIN, NX

      INTRINSIC          CHAR, ICHAR, INT, MIN, REAL

      INTEGER            IEEECK
      EXTERNAL           IEEECK

      GO TO ( 100, 100, 100, 400, 500, 600, 700, 800, 900, 1000,&
             1100 ) ISPEC

      ILAENV = -1
      RETURN

  100 CONTINUE

      ILAENV = 1
      SUBNAM = NAME
      IC = ICHAR( SUBNAM( 1:1 ) )
      IZ = ICHAR( 'Z' )
      IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN

         IF( IC.GE.97 .AND. IC.LE.122 ) THEN
            SUBNAM( 1:1 ) = CHAR( IC-32 )
            DO 10 I = 2, 6
               IC = ICHAR( SUBNAM( I:I ) )
               IF( IC.GE.97 .AND. IC.LE.122 )&
                 SUBNAM( I:I ) = CHAR( IC-32 )
   10       CONTINUE
         END IF

      ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN

         IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.&
            ( IC.GE.145 .AND. IC.LE.153 ) .OR.&
            ( IC.GE.162 .AND. IC.LE.169 ) ) THEN
            SUBNAM( 1:1 ) = CHAR( IC+64 )
            DO 20 I = 2, 6
               IC = ICHAR( SUBNAM( I:I ) )
               IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.&
                  ( IC.GE.145 .AND. IC.LE.153 ) .OR.&
                  ( IC.GE.162 .AND. IC.LE.169 ) )&
                 SUBNAM( I:I ) = CHAR( IC+64 )
   20       CONTINUE
         END IF

      ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN

         IF( IC.GE.225 .AND. IC.LE.250 ) THEN
            SUBNAM( 1:1 ) = CHAR( IC-32 )
            DO 30 I = 2, 6
               IC = ICHAR( SUBNAM( I:I ) )
               IF( IC.GE.225 .AND. IC.LE.250 )&
                 SUBNAM( I:I ) = CHAR( IC-32 )
   30       CONTINUE
         END IF
      END IF

      C1 = SUBNAM( 1:1 )
      SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
      CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
      IF( .NOT.( CNAME .OR. SNAME ) )&
        RETURN
      C2 = SUBNAM( 2:3 )
      C3 = SUBNAM( 4:6 )
      C4 = C3( 2:3 )

      GO TO ( 110, 200, 300 ) ISPEC

  110 CONTINUE

      NB = 1

      IF( C2.EQ.'GE' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.&
                 C3.EQ.'QLF' ) THEN
            IF( SNAME ) THEN
               NB = 32
            ELSE
               NB = 32
            END IF
         ELSE IF( C3.EQ.'HRD' ) THEN
            IF( SNAME ) THEN
               NB = 32
            ELSE
               NB = 32
            END IF
         ELSE IF( C3.EQ.'BRD' ) THEN
            IF( SNAME ) THEN
               NB = 32
            ELSE
               NB = 32
            END IF
         ELSE IF( C3.EQ.'TRI' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( C2.EQ.'PO' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( C2.EQ.'SY' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
            NB = 32
         ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN
            NB = 64
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            NB = 64
         ELSE IF( C3.EQ.'TRD' ) THEN
            NB = 32
         ELSE IF( C3.EQ.'GST' ) THEN
            NB = 64
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         END IF
      ELSE IF( C2.EQ.'GB' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               IF( N4.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            ELSE
               IF( N4.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            END IF
         END IF
      ELSE IF( C2.EQ.'PB' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               IF( N2.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            ELSE
               IF( N2.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            END IF
         END IF
      ELSE IF( C2.EQ.'TR' ) THEN
         IF( C3.EQ.'TRI' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( C2.EQ.'LA' ) THEN
         IF( C3.EQ.'UUM' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
         IF( C3.EQ.'EBZ' ) THEN
            NB = 1
         END IF
      END IF
      ILAENV = NB
      RETURN

  200 CONTINUE

      NBMIN = 2
      IF( C2.EQ.'GE' ) THEN
         IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.&
            C3.EQ.'QLF' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         ELSE IF( C3.EQ.'HRD' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         ELSE IF( C3.EQ.'BRD' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         ELSE IF( C3.EQ.'TRI' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         END IF
      ELSE IF( C2.EQ.'SY' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NBMIN = 8
            ELSE
               NBMIN = 8
            END IF
         ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
            NBMIN = 2
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
         IF( C3.EQ.'TRD' ) THEN
            NBMIN = 2
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         END IF
      END IF
      ILAENV = NBMIN
      RETURN

  300 CONTINUE

      NX = 0
      IF( C2.EQ.'GE' ) THEN
         IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.&
            C3.EQ.'QLF' ) THEN
            IF( SNAME ) THEN
               NX = 128
            ELSE
               NX = 128
            END IF
         ELSE IF( C3.EQ.'HRD' ) THEN
            IF( SNAME ) THEN
               NX = 128
            ELSE
               NX = 128
            END IF
         ELSE IF( C3.EQ.'BRD' ) THEN
            IF( SNAME ) THEN
               NX = 128
            ELSE
               NX = 128
            END IF
         END IF
      ELSE IF( C2.EQ.'SY' ) THEN
         IF( SNAME .AND. C3.EQ.'TRD' ) THEN
            NX = 32
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
         IF( C3.EQ.'TRD' ) THEN
            NX = 32
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NX = 128
            END IF
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.&
               C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.&
               C4.EQ.'BR' ) THEN
               NX = 128
            END IF
         END IF
      END IF
      ILAENV = NX
      RETURN

  400 CONTINUE

      ILAENV = 6
      RETURN

  500 CONTINUE

      ILAENV = 2