dlatms.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,040 行 · 第 1/3 页
F
1,040 行
ELSE
IPACK = -1
END IF
*
* Set certain internal parameters
*
MNMIN = MIN( M, N )
LLB = MIN( KL, M-1 )
UUB = MIN( KU, N-1 )
MR = MIN( M, N+LLB )
NC = MIN( N, M+UUB )
*
IF( IPACK.EQ.5 .OR. IPACK.EQ.6 ) THEN
MINLDA = UUB + 1
ELSE IF( IPACK.EQ.7 ) THEN
MINLDA = LLB + UUB + 1
ELSE
MINLDA = M
END IF
*
* Use Givens rotation method if bandwidth small enough,
* or if LDA is too small to store the matrix unpacked.
*
GIVENS = .FALSE.
IF( ISYM.EQ.1 ) THEN
IF( DBLE( LLB+UUB ).LT.0.3D0*DBLE( MAX( 1, MR+NC ) ) )
$ GIVENS = .TRUE.
ELSE
IF( 2*LLB.LT.M )
$ GIVENS = .TRUE.
END IF
IF( LDA.LT.M .AND. LDA.GE.MINLDA )
$ GIVENS = .TRUE.
*
* Set INFO if an error
*
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( M.NE.N .AND. ISYM.NE.1 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( IDIST.EQ.-1 ) THEN
INFO = -3
ELSE IF( ISYM.EQ.-1 ) THEN
INFO = -5
ELSE IF( ABS( MODE ).GT.6 ) THEN
INFO = -7
ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
$ THEN
INFO = -8
ELSE IF( KL.LT.0 ) THEN
INFO = -10
ELSE IF( KU.LT.0 .OR. ( ISYM.NE.1 .AND. KL.NE.KU ) ) THEN
INFO = -11
ELSE IF( IPACK.EQ.-1 .OR. ( ISYMPK.EQ.1 .AND. ISYM.EQ.1 ) .OR.
$ ( ISYMPK.EQ.2 .AND. ISYM.EQ.1 .AND. KL.GT.0 ) .OR.
$ ( ISYMPK.EQ.3 .AND. ISYM.EQ.1 .AND. KU.GT.0 ) .OR.
$ ( ISYMPK.NE.0 .AND. M.NE.N ) ) THEN
INFO = -12
ELSE IF( LDA.LT.MAX( 1, MINLDA ) ) THEN
INFO = -14
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DLATMS', -INFO )
RETURN
END IF
*
* Initialize random number generator
*
DO 10 I = 1, 4
ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
10 CONTINUE
*
IF( MOD( ISEED( 4 ), 2 ).NE.1 )
$ ISEED( 4 ) = ISEED( 4 ) + 1
*
* 2) Set up D if indicated.
*
* Compute D according to COND and MODE
*
CALL DLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = 1
RETURN
END IF
*
* Choose Top-Down if D is (apparently) increasing,
* Bottom-Up if D is (apparently) decreasing.
*
IF( ABS( D( 1 ) ).LE.ABS( D( MNMIN ) ) ) THEN
TOPDWN = .TRUE.
ELSE
TOPDWN = .FALSE.
END IF
*
IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
*
* Scale by DMAX
*
TEMP = ABS( D( 1 ) )
DO 20 I = 2, MNMIN
TEMP = MAX( TEMP, ABS( D( I ) ) )
20 CONTINUE
*
IF( TEMP.GT.ZERO ) THEN
ALPHA = DMAX / TEMP
ELSE
INFO = 2
RETURN
END IF
*
CALL DSCAL( MNMIN, ALPHA, D, 1 )
*
END IF
*
* 3) Generate Banded Matrix using Givens rotations.
* Also the special case of UUB=LLB=0
*
* Compute Addressing constants to cover all
* storage formats. Whether GE, SY, GB, or SB,
* upper or lower triangle or both,
* the (i,j)-th element is in
* A( i - ISKEW*j + IOFFST, j )
*
IF( IPACK.GT.4 ) THEN
ILDA = LDA - 1
ISKEW = 1
IF( IPACK.GT.5 ) THEN
IOFFST = UUB + 1
ELSE
IOFFST = 1
END IF
ELSE
ILDA = LDA
ISKEW = 0
IOFFST = 0
END IF
*
* IPACKG is the format that the matrix is generated in. If this is
* different from IPACK, then the matrix must be repacked at the
* end. It also signals how to compute the norm, for scaling.
*
IPACKG = 0
CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
*
* Diagonal Matrix -- We are done, unless it
* is to be stored SP/PP/TP (PACK='R' or 'C')
*
IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN
CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
IF( IPACK.LE.2 .OR. IPACK.GE.5 )
$ IPACKG = IPACK
*
ELSE IF( GIVENS ) THEN
*
* Check whether to use Givens rotations,
* Householder transformations, or nothing.
*
IF( ISYM.EQ.1 ) THEN
*
* Non-symmetric -- A = U D V
*
IF( IPACK.GT.4 ) THEN
IPACKG = IPACK
ELSE
IPACKG = 0
END IF
*
CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
*
IF( TOPDWN ) THEN
JKL = 0
DO 50 JKU = 1, UUB
*
* Transform from bandwidth JKL, JKU-1 to JKL, JKU
*
* Last row actually rotated is M
* Last column actually rotated is MIN( M+JKU, N )
*
DO 40 JR = 1, MIN( M+JKU, N ) + JKL - 1
EXTRA = ZERO
ANGLE = TWOPI*DLARND( 1, ISEED )
C = COS( ANGLE )
S = SIN( ANGLE )
ICOL = MAX( 1, JR-JKL )
IF( JR.LT.M ) THEN
IL = MIN( N, JR+JKU ) + 1 - ICOL
CALL DLAROT( .TRUE., JR.GT.JKL, .FALSE., IL, C,
$ S, A( JR-ISKEW*ICOL+IOFFST, ICOL ),
$ ILDA, EXTRA, DUMMY )
END IF
*
* Chase "EXTRA" back up
*
IR = JR
IC = ICOL
DO 30 JCH = JR - JKL, 1, -JKL - JKU
IF( IR.LT.M ) THEN
CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
$ IC+1 ), EXTRA, C, S, DUMMY )
END IF
IROW = MAX( 1, JCH-JKU )
IL = IR + 2 - IROW
TEMP = ZERO
ILTEMP = JCH.GT.JKU
CALL DLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S,
$ A( IROW-ISKEW*IC+IOFFST, IC ),
$ ILDA, TEMP, EXTRA )
IF( ILTEMP ) THEN
CALL DLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
$ IC+1 ), TEMP, C, S, DUMMY )
ICOL = MAX( 1, JCH-JKU-JKL )
IL = IC + 2 - ICOL
EXTRA = ZERO
CALL DLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE.,
$ IL, C, -S, A( IROW-ISKEW*ICOL+
$ IOFFST, ICOL ), ILDA, EXTRA,
$ TEMP )
IC = ICOL
IR = IROW
END IF
30 CONTINUE
40 CONTINUE
50 CONTINUE
*
JKU = UUB
DO 80 JKL = 1, LLB
*
* Transform from bandwidth JKL-1, JKU to JKL, JKU
*
DO 70 JC = 1, MIN( N+JKL, M ) + JKU - 1
EXTRA = ZERO
ANGLE = TWOPI*DLARND( 1, ISEED )
C = COS( ANGLE )
S = SIN( ANGLE )
IROW = MAX( 1, JC-JKU )
IF( JC.LT.N ) THEN
IL = MIN( M, JC+JKL ) + 1 - IROW
CALL DLAROT( .FALSE., JC.GT.JKU, .FALSE., IL, C,
$ S, A( IROW-ISKEW*JC+IOFFST, JC ),
$ ILDA, EXTRA, DUMMY )
END IF
*
* Chase "EXTRA" back up
*
IC = JC
IR = IROW
DO 60 JCH = JC - JKU, 1, -JKL - JKU
IF( IC.LT.N ) THEN
CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
$ IC+1 ), EXTRA, C, S, DUMMY )
END IF
ICOL = MAX( 1, JCH-JKL )
IL = IC + 2 - ICOL
TEMP = ZERO
ILTEMP = JCH.GT.JKL
CALL DLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S,
$ A( IR-ISKEW*ICOL+IOFFST, ICOL ),
$ ILDA, TEMP, EXTRA )
IF( ILTEMP ) THEN
CALL DLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
$ ICOL+1 ), TEMP, C, S, DUMMY )
IROW = MAX( 1, JCH-JKL-JKU )
IL = IR + 2 - IROW
EXTRA = ZERO
CALL DLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE.,
$ IL, C, -S, A( IROW-ISKEW*ICOL+
$ IOFFST, ICOL ), ILDA, EXTRA,
$ TEMP )
IC = ICOL
IR = IROW
END IF
60 CONTINUE
70 CONTINUE
80 CONTINUE
*
ELSE
*
* Bottom-Up -- Start at the bottom right.
*
JKL = 0
DO 110 JKU = 1, UUB
*
* Transform from bandwidth JKL, JKU-1 to JKL, JKU
*
* First row actually rotated is M
* First column actually rotated is MIN( M+JKU, N )
*
IENDCH = MIN( M, N+JKL ) - 1
DO 100 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1
EXTRA = ZERO
ANGLE = TWOPI*DLARND( 1, ISEED )
C = COS( ANGLE )
S = SIN( ANGLE )
IROW = MAX( 1, JC-JKU+1 )
IF( JC.GT.0 ) THEN
IL = MIN( M, JC+JKL+1 ) + 1 - IROW
CALL DLAROT( .FALSE., .FALSE., JC+JKL.LT.M, IL,
$ C, S, A( IROW-ISKEW*JC+IOFFST,
$ JC ), ILDA, DUMMY, EXTRA )
END IF
*
* Chase "EXTRA" back down
*
IC = JC
DO 90 JCH = JC + JKL, IENDCH, JKL + JKU
ILEXTR = IC.GT.0
IF( ILEXTR ) THEN
CALL DLARTG( A( JCH-ISKEW*IC+IOFFST, IC ),
$ EXTRA, C, S, DUMMY )
END IF
IC = MAX( 1, IC )
ICOL = MIN( N-1, JCH+JKU )
ILTEMP = JCH + JKU.LT.N
TEMP = ZERO
CALL DLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
$ C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
$ ILDA, EXTRA, TEMP )
IF( ILTEMP ) THEN
CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFST,
$ ICOL ), TEMP, C, S, DUMMY )
IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
EXTRA = ZERO
CALL DLAROT( .FALSE., .TRUE.,
$ JCH+JKL+JKU.LE.IENDCH, IL, C, S,
$ A( JCH-ISKEW*ICOL+IOFFST,
$ ICOL ), ILDA, TEMP, EXTRA )
IC = ICOL
END IF
90 CONTINUE
100 CONTINUE
110 CONTINUE
*
JKU = UUB
DO 140 JKL = 1, LLB
*
* Transform from bandwidth JKL-1, JKU to JKL, JKU
*
* First row actually rotated is MIN( N+JKL, M )
* First column actually rotated is N
*
IENDCH = MIN( N, M+JKU ) - 1
DO 130 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1
EXTRA = ZERO
ANGLE = TWOPI*DLARND( 1, ISEED )
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?