clatms.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,167 行 · 第 1/3 页
F
1,167 行
$ JCH+JKL+JKU.LE.IENDCH, IL, C, S,
$ A( IROW-ISKEW*JCH+IOFFST, JCH ),
$ ILDA, CTEMP, EXTRA )
IR = IROW
END IF
140 CONTINUE
150 CONTINUE
160 CONTINUE
*
END IF
*
ELSE
*
* Symmetric -- A = U D U'
* Hermitian -- A = U D U*
*
IPACKG = IPACK
IOFFG = IOFFST
*
IF( TOPDWN ) THEN
*
* Top-Down -- Generate Upper triangle only
*
IF( IPACK.GE.5 ) THEN
IPACKG = 6
IOFFG = UUB + 1
ELSE
IPACKG = 1
END IF
*
DO 170 J = 1, MNMIN
A( ( 1-ISKEW )*J+IOFFG, J ) = CMPLX( D( J ) )
170 CONTINUE
*
DO 200 K = 1, UUB
DO 190 JC = 1, N - 1
IROW = MAX( 1, JC-K )
IL = MIN( JC+1, K+2 )
EXTRA = CZERO
CTEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 )
ANGLE = TWOPI*SLARND( 1, ISEED )
C = COS( ANGLE )*CLARND( 5, ISEED )
S = SIN( ANGLE )*CLARND( 5, ISEED )
IF( CSYM ) THEN
CT = C
ST = S
ELSE
CTEMP = CONJG( CTEMP )
CT = CONJG( C )
ST = CONJG( S )
END IF
CALL CLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S,
$ A( IROW-ISKEW*JC+IOFFG, JC ), ILDA,
$ EXTRA, CTEMP )
CALL CLAROT( .TRUE., .TRUE., .FALSE.,
$ MIN( K, N-JC )+1, CT, ST,
$ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
$ CTEMP, DUMMY )
*
* Chase EXTRA back up the matrix
*
ICOL = JC
DO 180 JCH = JC - K, 1, -K
CALL CLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG,
$ ICOL+1 ), EXTRA, REALC, S, DUMMY )
DUMMY = CLARND( 5, ISEED )
C = CONJG( REALC*DUMMY )
S = CONJG( -S*DUMMY )
CTEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
IF( CSYM ) THEN
CT = C
ST = S
ELSE
CTEMP = CONJG( CTEMP )
CT = CONJG( C )
ST = CONJG( S )
END IF
CALL CLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
$ A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
$ ILDA, CTEMP, EXTRA )
IROW = MAX( 1, JCH-K )
IL = MIN( JCH+1, K+2 )
EXTRA = CZERO
CALL CLAROT( .FALSE., JCH.GT.K, .TRUE., IL, CT,
$ ST, A( IROW-ISKEW*JCH+IOFFG, JCH ),
$ ILDA, EXTRA, CTEMP )
ICOL = JCH
180 CONTINUE
190 CONTINUE
200 CONTINUE
*
* If we need lower triangle, copy from upper. Note that
* the order of copying is chosen to work for 'q' -> 'b'
*
IF( IPACK.NE.IPACKG .AND. IPACK.NE.3 ) THEN
DO 230 JC = 1, N
IROW = IOFFST - ISKEW*JC
IF( CSYM ) THEN
DO 210 JR = JC, MIN( N, JC+UUB )
A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
210 CONTINUE
ELSE
DO 220 JR = JC, MIN( N, JC+UUB )
A( JR+IROW, JC ) = CONJG( A( JC-ISKEW*JR+
$ IOFFG, JR ) )
220 CONTINUE
END IF
230 CONTINUE
IF( IPACK.EQ.5 ) THEN
DO 250 JC = N - UUB + 1, N
DO 240 JR = N + 2 - JC, UUB + 1
A( JR, JC ) = CZERO
240 CONTINUE
250 CONTINUE
END IF
IF( IPACKG.EQ.6 ) THEN
IPACKG = IPACK
ELSE
IPACKG = 0
END IF
END IF
ELSE
*
* Bottom-Up -- Generate Lower triangle only
*
IF( IPACK.GE.5 ) THEN
IPACKG = 5
IF( IPACK.EQ.6 )
$ IOFFG = 1
ELSE
IPACKG = 2
END IF
*
DO 260 J = 1, MNMIN
A( ( 1-ISKEW )*J+IOFFG, J ) = CMPLX( D( J ) )
260 CONTINUE
*
DO 290 K = 1, UUB
DO 280 JC = N - 1, 1, -1
IL = MIN( N+1-JC, K+2 )
EXTRA = CZERO
CTEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC )
ANGLE = TWOPI*SLARND( 1, ISEED )
C = COS( ANGLE )*CLARND( 5, ISEED )
S = SIN( ANGLE )*CLARND( 5, ISEED )
IF( CSYM ) THEN
CT = C
ST = S
ELSE
CTEMP = CONJG( CTEMP )
CT = CONJG( C )
ST = CONJG( S )
END IF
CALL CLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S,
$ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
$ CTEMP, EXTRA )
ICOL = MAX( 1, JC-K+1 )
CALL CLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL,
$ CT, ST, A( JC-ISKEW*ICOL+IOFFG,
$ ICOL ), ILDA, DUMMY, CTEMP )
*
* Chase EXTRA back down the matrix
*
ICOL = JC
DO 270 JCH = JC + K, N - 1, K
CALL CLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
$ EXTRA, REALC, S, DUMMY )
DUMMY = CLARND( 5, ISEED )
C = REALC*DUMMY
S = S*DUMMY
CTEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
IF( CSYM ) THEN
CT = C
ST = S
ELSE
CTEMP = CONJG( CTEMP )
CT = CONJG( C )
ST = CONJG( S )
END IF
CALL CLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
$ A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
$ ILDA, EXTRA, CTEMP )
IL = MIN( N+1-JCH, K+2 )
EXTRA = CZERO
CALL CLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL,
$ CT, ST, A( ( 1-ISKEW )*JCH+IOFFG,
$ JCH ), ILDA, CTEMP, EXTRA )
ICOL = JCH
270 CONTINUE
280 CONTINUE
290 CONTINUE
*
* If we need upper triangle, copy from lower. Note that
* the order of copying is chosen to work for 'b' -> 'q'
*
IF( IPACK.NE.IPACKG .AND. IPACK.NE.4 ) THEN
DO 320 JC = N, 1, -1
IROW = IOFFST - ISKEW*JC
IF( CSYM ) THEN
DO 300 JR = JC, MAX( 1, JC-UUB ), -1
A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
300 CONTINUE
ELSE
DO 310 JR = JC, MAX( 1, JC-UUB ), -1
A( JR+IROW, JC ) = CONJG( A( JC-ISKEW*JR+
$ IOFFG, JR ) )
310 CONTINUE
END IF
320 CONTINUE
IF( IPACK.EQ.6 ) THEN
DO 340 JC = 1, UUB
DO 330 JR = 1, UUB + 1 - JC
A( JR, JC ) = CZERO
330 CONTINUE
340 CONTINUE
END IF
IF( IPACKG.EQ.5 ) THEN
IPACKG = IPACK
ELSE
IPACKG = 0
END IF
END IF
END IF
*
* Ensure that the diagonal is real if Hermitian
*
IF( .NOT.CSYM ) THEN
DO 350 JC = 1, N
IROW = IOFFST + ( 1-ISKEW )*JC
A( IROW, JC ) = CMPLX( REAL( A( IROW, JC ) ) )
350 CONTINUE
END IF
*
END IF
*
ELSE
*
* 4) Generate Banded Matrix by first
* Rotating by random Unitary matrices,
* then reducing the bandwidth using Householder
* transformations.
*
* Note: we should get here only if LDA .ge. N
*
IF( ISYM.EQ.1 ) THEN
*
* Non-symmetric -- A = U D V
*
CALL CLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK,
$ IINFO )
ELSE
*
* Symmetric -- A = U D U' or
* Hermitian -- A = U D U*
*
IF( CSYM ) THEN
CALL CLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
ELSE
CALL CLAGHE( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
END IF
END IF
*
IF( IINFO.NE.0 ) THEN
INFO = 3
RETURN
END IF
END IF
*
* 5) Pack the matrix
*
IF( IPACK.NE.IPACKG ) THEN
IF( IPACK.EQ.1 ) THEN
*
* 'U' -- Upper triangular, not packed
*
DO 370 J = 1, M
DO 360 I = J + 1, M
A( I, J ) = CZERO
360 CONTINUE
370 CONTINUE
*
ELSE IF( IPACK.EQ.2 ) THEN
*
* 'L' -- Lower triangular, not packed
*
DO 390 J = 2, M
DO 380 I = 1, J - 1
A( I, J ) = CZERO
380 CONTINUE
390 CONTINUE
*
ELSE IF( IPACK.EQ.3 ) THEN
*
* 'C' -- Upper triangle packed Columnwise.
*
ICOL = 1
IROW = 0
DO 410 J = 1, M
DO 400 I = 1, J
IROW = IROW + 1
IF( IROW.GT.LDA ) THEN
IROW = 1
ICOL = ICOL + 1
END IF
A( IROW, ICOL ) = A( I, J )
400 CONTINUE
410 CONTINUE
*
ELSE IF( IPACK.EQ.4 ) THEN
*
* 'R' -- Lower triangle packed Columnwise.
*
ICOL = 1
IROW = 0
DO 430 J = 1, M
DO 420 I = J, M
IROW = IROW + 1
IF( IROW.GT.LDA ) THEN
IROW = 1
ICOL = ICOL + 1
END IF
A( IROW, ICOL ) = A( I, J )
420 CONTINUE
430 CONTINUE
*
ELSE IF( IPACK.GE.5 ) THEN
*
* 'B' -- The lower triangle is packed as a band matrix.
* 'Q' -- The upper triangle is packed as a band matrix.
* 'Z' -- The whole matrix is packed as a band matrix.
*
IF( IPACK.EQ.5 )
$ UUB = 0
IF( IPACK.EQ.6 )
$ LLB = 0
*
DO 450 J = 1, UUB
DO 440 I = MIN( J+LLB, M ), 1, -1
A( I-J+UUB+1, J ) = A( I, J )
440 CONTINUE
450 CONTINUE
*
DO 470 J = UUB + 2, N
DO 460 I = J - UUB, MIN( J+LLB, M )
A( I-J+UUB+1, J ) = A( I, J )
460 CONTINUE
470 CONTINUE
END IF
*
* If packed, zero out extraneous elements.
*
* Symmetric/Triangular Packed --
* zero out everything after A(IROW,ICOL)
*
IF( IPACK.EQ.3 .OR. IPACK.EQ.4 ) THEN
DO 490 JC = ICOL, M
DO 480 JR = IROW + 1, LDA
A( JR, JC ) = CZERO
480 CONTINUE
IROW = 0
490 CONTINUE
*
ELSE IF( IPACK.GE.5 ) THEN
*
* Packed Band --
* 1st row is now in A( UUB+2-j, j), zero above it
* m-th row is now in A( M+UUB-j,j), zero below it
* last non-zero diagonal is now in A( UUB+LLB+1,j ),
* zero below it, too.
*
IR1 = UUB + LLB + 2
IR2 = UUB + M + 2
DO 520 JC = 1, N
DO 500 JR = 1, UUB + 1 - JC
A( JR, JC ) = CZERO
500 CONTINUE
DO 510 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
A( JR, JC ) = CZERO
510 CONTINUE
520 CONTINUE
END IF
END IF
*
RETURN
*
* End of CLATMS
*
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?