matrix.f90

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

      RETURN

  600 CONTINUE 

      ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
      RETURN

  700 CONTINUE

      ILAENV = 1
      RETURN

  800 CONTINUE

      ILAENV = 50
      RETURN

  900 CONTINUE

      ILAENV = 25
      RETURN

 1000 CONTINUE

      ILAENV = 0
      RETURN

 1100 CONTINUE

      ILAENV = 0
      RETURN

      END
!==========================================================================

      LOGICAL          FUNCTION LSAME( CA, CB )

      CHARACTER          CA, CB


      INTRINSIC          ICHAR

      INTEGER            INTA, INTB, ZCODE

      LSAME = CA.EQ.CB
      IF( LSAME )&
        RETURN

      ZCODE = ICHAR( 'Z' )

      INTA = ICHAR( CA )
      INTB = ICHAR( CB )

      IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN

         IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
         IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32

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

         IF( INTA.GE.129 .AND. INTA.LE.137 .OR.&
            INTA.GE.145 .AND. INTA.LE.153 .OR.&
            INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64
         IF( INTB.GE.129 .AND. INTB.LE.137 .OR.&
            INTB.GE.145 .AND. INTB.LE.153 .OR.&
            INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64

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

         IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
         IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
      END IF
      LSAME = INTA.EQ.INTB

      END
!==========================================================================


      subroutine saxpy(n,sa,sx,incx,sy,incy)

      real sx(*),sy(*),sa
      integer i,incx,incy,ix,iy,m,mp1,n

      if(n.le.0)return
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20

      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        sy(iy) = sy(iy) + sa*sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return

   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sy(i) = sy(i) + sa*sx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        sy(i) = sy(i) + sa*sx(i)
        sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
        sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
        sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
   50 continue
      return
      end
!========================================================================

      real function sdot(n,sx,incx,sy,incy)

      real sx(*),sy(*),stemp
      integer i,incx,incy,ix,iy,m,mp1,n

      stemp = 0.0e0
      sdot = 0.0e0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20

      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        stemp = stemp + sx(ix)*sy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      sdot = stemp
      return

   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        stemp = stemp + sx(i)*sy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +&
        sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
   50 continue
   60 sdot = stemp
      return
      end
!===========================================================================
      SUBROUTINE SGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,&
                        BETA, Y, INCY )

      REAL               ALPHA, BETA
      INTEGER            INCX, INCY, LDA, M, N
      CHARACTER*1        TRANS

      REAL               A( LDA, * ), X( * ), Y( * )




      REAL               ONE         , ZERO
      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )

      REAL               TEMP
      INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY

      LOGICAL            LSAME
      EXTERNAL           LSAME

      EXTERNAL           XERBLA

      INTRINSIC          MAX

      INFO = 0
      IF     ( .NOT.LSAME( TRANS, 'N' ).AND.&
              .NOT.LSAME( TRANS, 'T' ).AND.&
              .NOT.LSAME( TRANS, 'C' )      )THEN
         INFO = 1
      ELSE IF( M.LT.0 )THEN
         INFO = 2
      ELSE IF( N.LT.0 )THEN
         INFO = 3
      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
         INFO = 6
      ELSE IF( INCX.EQ.0 )THEN
         INFO = 8
      ELSE IF( INCY.EQ.0 )THEN
         INFO = 11
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'SGEMV ', INFO )
         RETURN
      END IF

      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.&
         ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )&
        RETURN

      IF( LSAME( TRANS, 'N' ) )THEN
         LENX = N
         LENY = M
      ELSE
         LENX = M
         LENY = N
      END IF
      IF( INCX.GT.0 )THEN
         KX = 1
      ELSE
         KX = 1 - ( LENX - 1 )*INCX
      END IF
      IF( INCY.GT.0 )THEN
         KY = 1
      ELSE
         KY = 1 - ( LENY - 1 )*INCY
      END IF

      IF( BETA.NE.ONE )THEN
         IF( INCY.EQ.1 )THEN
            IF( BETA.EQ.ZERO )THEN
               DO 10, I = 1, LENY
                  Y( I ) = ZERO
   10          CONTINUE
            ELSE
               DO 20, I = 1, LENY
                  Y( I ) = BETA*Y( I )
   20          CONTINUE
            END IF
         ELSE
            IY = KY
            IF( BETA.EQ.ZERO )THEN
               DO 30, I = 1, LENY
                  Y( IY ) = ZERO
                  IY      = IY   + INCY
   30          CONTINUE
            ELSE
               DO 40, I = 1, LENY
                  Y( IY ) = BETA*Y( IY )
                  IY      = IY           + INCY
   40          CONTINUE
            END IF
         END IF
      END IF
      IF( ALPHA.EQ.ZERO )&
        RETURN
      IF( LSAME( TRANS, 'N' ) )THEN

         JX = KX
         IF( INCY.EQ.1 )THEN
            DO 60, J = 1, N
               IF( X( JX ).NE.ZERO )THEN
                  TEMP = ALPHA*X( JX )
                  DO 50, I = 1, M
                     Y( I ) = Y( I ) + TEMP*A( I, J )
   50             CONTINUE
               END IF
               JX = JX + INCX
   60       CONTINUE
         ELSE
            DO 80, J = 1, N
               IF( X( JX ).NE.ZERO )THEN
                  TEMP = ALPHA*X( JX )
                  IY   = KY
                  DO 70, I = 1, M
                     Y( IY ) = Y( IY ) + TEMP*A( I, J )
                     IY      = IY      + INCY
   70             CONTINUE
               END IF
               JX = JX + INCX
   80       CONTINUE
         END IF
      ELSE

         JY = KY
         IF( INCX.EQ.1 )THEN
            DO 100, J = 1, N
               TEMP = ZERO
               DO 90, I = 1, M
                  TEMP = TEMP + A( I, J )*X( I )
   90          CONTINUE
               Y( JY ) = Y( JY ) + ALPHA*TEMP
               JY      = JY      + INCY
  100       CONTINUE
         ELSE
            DO 120, J = 1, N
               TEMP = ZERO
               IX   = KX
               DO 110, I = 1, M
                  TEMP = TEMP + A( I, J )*X( IX )
                  IX   = IX   + INCX
  110          CONTINUE
               Y( JY ) = Y( JY ) + ALPHA*TEMP
               JY      = JY      + INCY
  120       CONTINUE
         END IF
      END IF

      RETURN

      END
!==========================================================================

      SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,&
                        RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,&
                        NAB, WORK, IWORK, INFO )

      INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
      REAL               ABSTOL, PIVMIN, RELTOL

      INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
      REAL               AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),&
                        WORK( * )


      REAL               ZERO, TWO, HALF
      PARAMETER          ( ZERO = 0.0E0, TWO = 2.0E0,&
                        HALF = 1.0E0 / TWO )

      INTEGER            ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,&
                        KLNEW
      REAL               TMP1, TMP2

      INTRINSIC          ABS, MAX, MIN

      INFO = 0
      IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
         INFO = -1
         RETURN
      END IF

      IF( IJOB.EQ.1 ) THEN

         MOUT = 0

         DO 30 JI = 1, MINP
            DO 20 JP = 1, 2
               TMP1 = D( 1 ) - AB( JI, JP )
               IF( ABS( TMP1 ).LT.PIVMIN )&
                 TMP1 = -PIVMIN
               NAB( JI, JP ) = 0
               IF( TMP1.LE.ZERO )&
                 NAB( JI, JP ) = 1

               DO 10 J = 2, N
                  TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
                  IF( ABS( TMP1 ).LT.PIVMIN )&
                    TMP1 = -PIVMIN
                  IF( TMP1.LE.ZERO )&
                    NAB( JI, JP ) = NAB( JI, JP ) + 1
   10          CONTINUE
   20       CONTINUE
            MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
   30    CONTINUE
         RETURN
      END IF

      KF = 1
      KL = MINP

      IF( IJOB.EQ.2 ) THEN
         DO 40 JI = 1, MINP
            C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
   40    CONTINUE
      END IF

      DO 130 JIT = 1, NITMAX

         IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN

            DO 60 JI = KF, KL

               WORK( JI ) = D( 1 ) - C( JI )
               IWORK( JI ) = 0
               IF( WORK( JI ).LE.PIVMIN ) THEN
                  IWORK( JI ) = 1
                  WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
               END IF

               DO 50 J = 2, N
                  WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
                  IF( WORK( JI ).LE.PIVMIN ) THEN
                     IWORK( JI ) = IWORK( JI ) + 1
                     WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
                  END IF
   50          CONTINUE
   60       CONTINUE

            IF( IJOB.LE.2 ) THEN

               KLNEW = KL
               DO 70 JI = KF, KL

                  IWORK( JI ) = MIN( NAB( JI, 2 ),&
                               MAX( NAB( JI, 1 ), IWORK( JI ) ) )


                  IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN

                     AB( JI, 2 ) = C( JI )

                  ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN

                     AB( JI, 1 ) = C( JI )
                  ELSE
                     KLNEW = KLNEW + 1
                     IF( KLNEW.LE.MMAX ) THEN


                        AB( KLNEW, 2 ) = AB( JI, 2 )
                        NAB( KLNEW, 2 ) = NAB( JI, 2 )
                        AB( KLNEW, 1 ) = C( JI )
                        NAB( KLNEW, 1 ) = IWORK( JI )
                        AB( JI, 2 ) = C( JI )
                        NAB( JI, 2 ) = IWORK( JI )
                     ELSE
                        INFO = MMAX + 1
                     END IF
                  END IF
   70          CONTINUE
               IF( INFO.NE.0 )&
                 RETURN
               KL = KLNEW
            ELSE

               DO 80 JI = KF, KL
                  IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
                     AB( JI, 1 ) = C( JI )
                     NAB( JI, 1 ) = IWORK( JI )
                  END IF
                  IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
                     AB( JI, 2 ) = C( JI )
                     NAB( JI, 2 ) = IWORK( JI )
                  END IF
   80          CONTINUE
            END IF

         ELSE

            KLNEW = KL
            DO 100 JI = KF, KL

               TMP1 = C( JI )
               TMP2 = D( 1 ) - TMP1
               ITMP1 = 0
               IF( TMP2.LE.PIVMIN ) THEN
                  ITMP1 = 1
                  TMP2 = MIN( TMP2, -PIVMIN )
               END IF

               DO 90 J = 2, N
                  TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
                  IF( TMP2.LE.PIVMIN ) THEN
                     ITMP1 = ITMP1 + 1
                     TMP2 = MIN( TMP2, -PIVMIN )
                  END IF
   90          CONTINUE

               IF( IJOB.LE.2 ) THEN

                  ITMP1 = MIN( NAB( JI, 2 ),&
                         MAX( NAB( JI, 1 ), ITMP1 ) )

                  IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN

                     AB( JI, 2 ) = TMP1

                  ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN

                     AB( JI, 1 ) = TMP1
                  ELSE IF( KLNEW.LT.MMAX ) THEN

                     KLNEW = KLNEW + 1
                     AB( KLNEW, 2 ) = AB( JI, 2 )
                     NAB( KLNEW, 2 ) = NAB( JI, 2 )
                     AB( KLNEW, 1 ) = TMP1
                     NAB( KLNEW, 1 ) = ITMP1
                     AB( JI, 2 ) = TMP1
                     NAB( JI, 2 ) = ITMP1
                  ELSE
                     INFO = MMAX + 1
                     RETURN
                  END IF
               ELSE

                  IF( ITMP1.LE.NVAL( JI ) ) THEN
                     AB( JI, 1 ) = TMP1
                     NAB( JI, 1 ) = ITMP1
                  END IF