clattb.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 617 行 · 第 1/2 页

*
*        Type 11:  Make the first diagonal element in the solve small to
*        cause immediate overflow when dividing by T(j,j).
*        In type 11, the offdiagonal elements are small (CNORM(j) < 1).
*
         CALL CLARNV( 2, ISEED, N, B )
         TSCAL = ONE / REAL( KD+1 )
         IF( UPPER ) THEN
            DO 140 J = 1, N
               LENJ = MIN( J-1, KD )
               IF( LENJ.GT.0 ) THEN
                  CALL CLARNV( 4, ISEED, LENJ, AB( KD+2-LENJ, J ) )
                  CALL CSSCAL( LENJ, TSCAL, AB( KD+2-LENJ, J ), 1 )
               END IF
               AB( KD+1, J ) = CLARND( 5, ISEED )
  140       CONTINUE
            AB( KD+1, N ) = SMLNUM*AB( KD+1, N )
         ELSE
            DO 150 J = 1, N
               LENJ = MIN( N-J, KD )
               IF( LENJ.GT.0 ) THEN
                  CALL CLARNV( 4, ISEED, LENJ, AB( 2, J ) )
                  CALL CSSCAL( LENJ, TSCAL, AB( 2, J ), 1 )
               END IF
               AB( 1, J ) = CLARND( 5, ISEED )
  150       CONTINUE
            AB( 1, 1 ) = SMLNUM*AB( 1, 1 )
         END IF
*
      ELSE IF( IMAT.EQ.12 ) THEN
*
*        Type 12:  Make the first diagonal element in the solve small to
*        cause immediate overflow when dividing by T(j,j).
*        In type 12, the offdiagonal elements are O(1) (CNORM(j) > 1).
*
         CALL CLARNV( 2, ISEED, N, B )
         IF( UPPER ) THEN
            DO 160 J = 1, N
               LENJ = MIN( J-1, KD )
               IF( LENJ.GT.0 )
     $            CALL CLARNV( 4, ISEED, LENJ, AB( KD+2-LENJ, J ) )
               AB( KD+1, J ) = CLARND( 5, ISEED )
  160       CONTINUE
            AB( KD+1, N ) = SMLNUM*AB( KD+1, N )
         ELSE
            DO 170 J = 1, N
               LENJ = MIN( N-J, KD )
               IF( LENJ.GT.0 )
     $            CALL CLARNV( 4, ISEED, LENJ, AB( 2, J ) )
               AB( 1, J ) = CLARND( 5, ISEED )
  170       CONTINUE
            AB( 1, 1 ) = SMLNUM*AB( 1, 1 )
         END IF
*
      ELSE IF( IMAT.EQ.13 ) THEN
*
*        Type 13:  T is diagonal with small numbers on the diagonal to
*        make the growth factor underflow, but a small right hand side
*        chosen so that the solution does not overflow.
*
         IF( UPPER ) THEN
            JCOUNT = 1
            DO 190 J = N, 1, -1
               DO 180 I = MAX( 1, KD+1-( J-1 ) ), KD
                  AB( I, J ) = ZERO
  180          CONTINUE
               IF( JCOUNT.LE.2 ) THEN
                  AB( KD+1, J ) = SMLNUM*CLARND( 5, ISEED )
               ELSE
                  AB( KD+1, J ) = CLARND( 5, ISEED )
               END IF
               JCOUNT = JCOUNT + 1
               IF( JCOUNT.GT.4 )
     $            JCOUNT = 1
  190       CONTINUE
         ELSE
            JCOUNT = 1
            DO 210 J = 1, N
               DO 200 I = 2, MIN( N-J+1, KD+1 )
                  AB( I, J ) = ZERO
  200          CONTINUE
               IF( JCOUNT.LE.2 ) THEN
                  AB( 1, J ) = SMLNUM*CLARND( 5, ISEED )
               ELSE
                  AB( 1, J ) = CLARND( 5, ISEED )
               END IF
               JCOUNT = JCOUNT + 1
               IF( JCOUNT.GT.4 )
     $            JCOUNT = 1
  210       CONTINUE
         END IF
*
*        Set the right hand side alternately zero and small.
*
         IF( UPPER ) THEN
            B( 1 ) = ZERO
            DO 220 I = N, 2, -2
               B( I ) = ZERO
               B( I-1 ) = SMLNUM*CLARND( 5, ISEED )
  220       CONTINUE
         ELSE
            B( N ) = ZERO
            DO 230 I = 1, N - 1, 2
               B( I ) = ZERO
               B( I+1 ) = SMLNUM*CLARND( 5, ISEED )
  230       CONTINUE
         END IF
*
      ELSE IF( IMAT.EQ.14 ) THEN
*
*        Type 14:  Make the diagonal elements small to cause gradual
*        overflow when dividing by T(j,j).  To control the amount of
*        scaling needed, the matrix is bidiagonal.
*
         TEXP = ONE / REAL( KD+1 )
         TSCAL = SMLNUM**TEXP
         CALL CLARNV( 4, ISEED, N, B )
         IF( UPPER ) THEN
            DO 250 J = 1, N
               DO 240 I = MAX( 1, KD+2-J ), KD
                  AB( I, J ) = ZERO
  240          CONTINUE
               IF( J.GT.1 .AND. KD.GT.0 )
     $            AB( KD, J ) = CMPLX( -ONE, -ONE )
               AB( KD+1, J ) = TSCAL*CLARND( 5, ISEED )
  250       CONTINUE
            B( N ) = CMPLX( ONE, ONE )
         ELSE
            DO 270 J = 1, N
               DO 260 I = 3, MIN( N-J+1, KD+1 )
                  AB( I, J ) = ZERO
  260          CONTINUE
               IF( J.LT.N .AND. KD.GT.0 )
     $            AB( 2, J ) = CMPLX( -ONE, -ONE )
               AB( 1, J ) = TSCAL*CLARND( 5, ISEED )
  270       CONTINUE
            B( 1 ) = CMPLX( ONE, ONE )
         END IF
*
      ELSE IF( IMAT.EQ.15 ) THEN
*
*        Type 15:  One zero diagonal element.
*
         IY = N / 2 + 1
         IF( UPPER ) THEN
            DO 280 J = 1, N
               LENJ = MIN( J, KD+1 )
               CALL CLARNV( 4, ISEED, LENJ, AB( KD+2-LENJ, J ) )
               IF( J.NE.IY ) THEN
                  AB( KD+1, J ) = CLARND( 5, ISEED )*TWO
               ELSE
                  AB( KD+1, J ) = ZERO
               END IF
  280       CONTINUE
         ELSE
            DO 290 J = 1, N
               LENJ = MIN( N-J+1, KD+1 )
               CALL CLARNV( 4, ISEED, LENJ, AB( 1, J ) )
               IF( J.NE.IY ) THEN
                  AB( 1, J ) = CLARND( 5, ISEED )*TWO
               ELSE
                  AB( 1, J ) = ZERO
               END IF
  290       CONTINUE
         END IF
         CALL CLARNV( 2, ISEED, N, B )
         CALL CSSCAL( N, TWO, B, 1 )
*
      ELSE IF( IMAT.EQ.16 ) THEN
*
*        Type 16:  Make the offdiagonal elements large to cause overflow
*        when adding a column of T.  In the non-transposed case, the
*        matrix is constructed to cause overflow when adding a column in
*        every other step.
*
         TSCAL = UNFL / ULP
         TSCAL = ( ONE-ULP ) / TSCAL
         DO 310 J = 1, N
            DO 300 I = 1, KD + 1
               AB( I, J ) = ZERO
  300       CONTINUE
  310    CONTINUE
         TEXP = ONE
         IF( KD.GT.0 ) THEN
            IF( UPPER ) THEN
               DO 330 J = N, 1, -KD
                  DO 320 I = J, MAX( 1, J-KD+1 ), -2
                     AB( 1+( J-I ), I ) = -TSCAL / REAL( KD+2 )
                     AB( KD+1, I ) = ONE
                     B( I ) = TEXP*( ONE-ULP )
                     IF( I.GT.MAX( 1, J-KD+1 ) ) THEN
                        AB( 2+( J-I ), I-1 ) = -( TSCAL / REAL( KD+2 ) )
     $                                          / REAL( KD+3 )
                        AB( KD+1, I-1 ) = ONE
                        B( I-1 ) = TEXP*REAL( ( KD+1 )*( KD+1 )+KD )
                     END IF
                     TEXP = TEXP*TWO
  320             CONTINUE
                  B( MAX( 1, J-KD+1 ) ) = ( REAL( KD+2 ) /
     $                                    REAL( KD+3 ) )*TSCAL
  330          CONTINUE
            ELSE
               DO 350 J = 1, N, KD
                  TEXP = ONE
                  LENJ = MIN( KD+1, N-J+1 )
                  DO 340 I = J, MIN( N, J+KD-1 ), 2
                     AB( LENJ-( I-J ), J ) = -TSCAL / REAL( KD+2 )
                     AB( 1, J ) = ONE
                     B( J ) = TEXP*( ONE-ULP )
                     IF( I.LT.MIN( N, J+KD-1 ) ) THEN
                        AB( LENJ-( I-J+1 ), I+1 ) = -( TSCAL /
     $                     REAL( KD+2 ) ) / REAL( KD+3 )
                        AB( 1, I+1 ) = ONE
                        B( I+1 ) = TEXP*REAL( ( KD+1 )*( KD+1 )+KD )
                     END IF
                     TEXP = TEXP*TWO
  340             CONTINUE
                  B( MIN( N, J+KD-1 ) ) = ( REAL( KD+2 ) /
     $                                    REAL( KD+3 ) )*TSCAL
  350          CONTINUE
            END IF
         END IF
*
      ELSE IF( IMAT.EQ.17 ) THEN
*
*        Type 17:  Generate a unit triangular matrix with elements
*        between -1 and 1, and make the right hand side large so that it
*        requires scaling.
*
         IF( UPPER ) THEN
            DO 360 J = 1, N
               LENJ = MIN( J-1, KD )
               CALL CLARNV( 4, ISEED, LENJ, AB( KD+1-LENJ, J ) )
               AB( KD+1, J ) = REAL( J )
  360       CONTINUE
         ELSE
            DO 370 J = 1, N
               LENJ = MIN( N-J, KD )
               IF( LENJ.GT.0 )
     $            CALL CLARNV( 4, ISEED, LENJ, AB( 2, J ) )
               AB( 1, J ) = REAL( J )
  370       CONTINUE
         END IF
*
*        Set the right hand side so that the largest value is BIGNUM.
*
         CALL CLARNV( 2, ISEED, N, B )
         IY = ICAMAX( N, B, 1 )
         BNORM = ABS( B( IY ) )
         BSCAL = BIGNUM / MAX( ONE, BNORM )
         CALL CSSCAL( N, BSCAL, B, 1 )
*
      ELSE IF( IMAT.EQ.18 ) THEN
*
*        Type 18:  Generate a triangular matrix with elements between
*        BIGNUM/(KD+1) and BIGNUM so that at least one of the column
*        norms will exceed BIGNUM.
*        1/3/91:  CLATBS no longer can handle this case
*
         TLEFT = BIGNUM / REAL( KD+1 )
         TSCAL = BIGNUM*( REAL( KD+1 ) / REAL( KD+2 ) )
         IF( UPPER ) THEN
            DO 390 J = 1, N
               LENJ = MIN( J, KD+1 )
               CALL CLARNV( 5, ISEED, LENJ, AB( KD+2-LENJ, J ) )
               CALL SLARNV( 1, ISEED, LENJ, RWORK( KD+2-LENJ ) )
               DO 380 I = KD + 2 - LENJ, KD + 1
                  AB( I, J ) = AB( I, J )*( TLEFT+RWORK( I )*TSCAL )
  380          CONTINUE
  390       CONTINUE
         ELSE
            DO 410 J = 1, N
               LENJ = MIN( N-J+1, KD+1 )
               CALL CLARNV( 5, ISEED, LENJ, AB( 1, J ) )
               CALL SLARNV( 1, ISEED, LENJ, RWORK )
               DO 400 I = 1, LENJ
                  AB( I, J ) = AB( I, J )*( TLEFT+RWORK( I )*TSCAL )
  400          CONTINUE
  410       CONTINUE
         END IF
         CALL CLARNV( 2, ISEED, N, B )
         CALL CSSCAL( N, TWO, B, 1 )
      END IF
*
*     Flip the matrix if the transpose will be used.
*
      IF( .NOT.LSAME( TRANS, 'N' ) ) THEN
         IF( UPPER ) THEN
            DO 420 J = 1, N / 2
               LENJ = MIN( N-2*J+1, KD+1 )
               CALL CSWAP( LENJ, AB( KD+1, J ), LDAB-1,
     $                     AB( KD+2-LENJ, N-J+1 ), -1 )
  420       CONTINUE
         ELSE
            DO 430 J = 1, N / 2
               LENJ = MIN( N-2*J+1, KD+1 )
               CALL CSWAP( LENJ, AB( 1, J ), 1, AB( LENJ, N-J+2-LENJ ),
     $                     -LDAB+1 )
  430       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     End of CLATTB
*
      END