!************************************************************************** !* LINEAR ALGEBRA LAPACK LIBRARY IN FORTRAN 90 (EXTRACT 3) * !* ---------------------------------------------------------------------- * !* Ref.: From LAPACK Fortran 77 Library, University of Tennessee, Univ. * !* of California Berkeley, NAG Ltd., Courant Institute, Argonne * !* National Lab, and Rice University, March 31, 1993 [BIBLI 17]. * !* * !* Fortran 90 Release 1.0 By J.-P. Moreau, Paris. * !************************************************************************** MODULE LAPACK3 !subroutines to solve a complex band linear system. CONTAINS SUBROUTINE CGBSV( N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO ) USE LAPACK2 ! ! -- LAPACK driver routine (version 3.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! .. Scalar Arguments .. INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS ! .. Array Arguments .. INTEGER IPIV( * ) COMPLEX AB( LDAB, * ), B( LDB, * ) ! .. ! Purpose ! ======= ! CGBSV computes the solution to a complex system of linear equations ! A * X = B, where A is a band matrix of order N with KL subdiagonals ! and KU superdiagonals, and X and B are N-by-NRHS matrices. ! The LU decomposition with partial pivoting and row interchanges is ! used to factor A as A = L * U, where L is a product of permutation ! and unit lower triangular matrices with KL subdiagonals, and U is ! upper triangular with KL+KU superdiagonals. The factored form of A ! is then used to solve the system of equations A * X = B. ! Arguments ! ========= ! N (input) INTEGER ! The number of linear equations, i.e., the order of the ! matrix A. N >= 0. ! KL (input) INTEGER ! The number of subdiagonals within the band of A. KL >= 0. ! KU (input) INTEGER ! The number of superdiagonals within the band of A. KU >= 0. ! NRHS (input) INTEGER ! The number of right hand sides, i.e., the number of columns ! of the matrix B. NRHS >= 0. ! AB (input/output) COMPLEX array, dimension (LDAB,N) ! On entry, the matrix A in band storage, in rows KL+1 to ! 2*KL+KU+1; rows 1 to KL of the array need not be set. ! The j-th column of A is stored in the j-th column of the ! array AB as follows: ! AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+KL) ! On exit, details of the factorization: U is stored as an ! upper triangular band matrix with KL+KU superdiagonals in ! rows 1 to KL+KU+1, and the multipliers used during the ! factorization are stored in rows KL+KU+2 to 2*KL+KU+1. ! See below for further details. ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= 2*KL+KU+1. ! IPIV (output) INTEGER array, dimension (N) ! The pivot indices that define the permutation matrix P; ! row i of the matrix was interchanged with row IPIV(i). ! B (input/output) COMPLEX array, dimension (LDB,NRHS) ! On entry, the N-by-NRHS right hand side matrix B. ! On exit, if INFO = 0, the N-by-NRHS solution matrix X. ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(1,N). ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! > 0: if INFO = i, U(i,i) is exactly zero. The factorization ! has been completed, but the factor U is exactly ! singular, and the solution has not been computed. ! Further Details ! =============== ! The band storage scheme is illustrated by the following example, when ! M = N = 6, KL = 2, KU = 1: ! On entry: On exit: ! * * * + + + * * * u14 u25 u36 ! * * + + + + * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * ! a31 a42 a53 a64 * * m31 m42 m53 m64 * * ! Array elements marked * are not used by the routine; elements marked ! + need not be set on entry, but are required by the routine to store ! elements of U because of fill-in resulting from the row interchanges. ! ====================================================================== ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. !* Test the input parameters. INFO = 0 IF( N.LT.0 ) THEN INFO = -1 ELSE IF( KL.LT.0 ) THEN INFO = -2 ELSE IF( KU.LT.0 ) THEN INFO = -3 ELSE IF( NRHS.LT.0 ) THEN INFO = -4 ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN INFO = -6 ELSE IF( LDB.LT.MAX( N, 1 ) ) THEN INFO = -9 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGBSV ', -INFO ) RETURN END IF ! Compute the LU factorization of the band matrix A. CALL CGBTRF( N, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( INFO.EQ.0 ) THEN ! Solve the system A*X = B, overwriting B with X. ! CALL CGBTRS( 'No transpose', N, KL, KU, NRHS, AB, LDAB, IPIV, & B, LDB, INFO ) END IF RETURN END SUBROUTINE CGBSV ! ----------------------------------------------------------------------- SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) USE LAPACK2 USE LAPACK1 ! -- LAPACK routine (version 3.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! September 30, 1994 ! .. Scalar Arguments .. INTEGER INFO, KL, KU, LDAB, M, N ! .. ! .. Array Arguments .. INTEGER IPIV( * ) COMPLEX AB( LDAB, * ) ! Purpose ! ======= ! CGBTRF computes an LU factorization of a complex m-by-n band matrix A ! using partial pivoting with row interchanges. ! This is the blocked version of the algorithm, calling Level 3 BLAS. ! Arguments ! ========= ! M (input) INTEGER ! The number of rows of the matrix A. M >= 0. ! N (input) INTEGER ! The number of columns of the matrix A. N >= 0. ! KL (input) INTEGER ! The number of subdiagonals within the band of A. KL >= 0. ! KU (input) INTEGER ! The number of superdiagonals within the band of A. KU >= 0. ! AB (input/output) COMPLEX array, dimension (LDAB,N) ! On entry, the matrix A in band storage, in rows KL+1 to ! 2*KL+KU+1; rows 1 to KL of the array need not be set. ! The j-th column of A is stored in the j-th column of the ! array AB as follows: ! AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) ! On exit, details of the factorization: U is stored as an ! upper triangular band matrix with KL+KU superdiagonals in ! rows 1 to KL+KU+1, and the multipliers used during the ! factorization are stored in rows KL+KU+2 to 2*KL+KU+1. ! See below for further details. ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= 2*KL+KU+1. ! IPIV (output) INTEGER array, dimension (min(M,N)) ! The pivot indices; for 1 <= i <= min(M,N), row i of the ! matrix was interchanged with row IPIV(i). ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! > 0: if INFO = +i, U(i,i) is exactly zero. The factorization ! has been completed, but the factor U is exactly ! singular, and division by zero will occur if it is used ! to solve a system of equations. ! Further Details ! =============== ! The band storage scheme is illustrated by the following example, when ! M = N = 6, KL = 2, KU = 1: ! On entry: On exit: ! * * * + + + * * * u14 u25 u36 ! * * + + + + * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * ! a31 a42 a53 a64 * * m31 m42 m53 m64 * * ! Array elements marked * are not used by the routine; elements marked ! + need not be set on entry, but are required by the routine to store ! elements of U because of fill-in resulting from the row interchanges. ! ====================================================================== ! .. Parameters .. COMPLEX ONE, ZERO PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), & ZERO = ( 0.0E+0, 0.0E+0 ) ) INTEGER NBMAX, LDWORK PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 ) ! .. Local Scalars .. INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, & JU, K2, KM, KV, NB, NW COMPLEX TEMP ! .. ! .. Local Arrays .. COMPLEX WORK13( LDWORK, NBMAX ), & WORK31( LDWORK, NBMAX ) ! .. ! .. External Functions .. ! INTEGER ICAMAX, ILAENV ! EXTERNAL ICAMAX, ILAENV ! .. ! .. External Subroutines .. ! EXTERNAL CCOPY, CGBTF2, CGEMM, CGERU, CLASWP, CSCAL, & ! CSWAP, CTRSM, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. Executable Statements .. ! KV is the number of superdiagonals in the factor U, allowing for ! fill-in KV = KU + KL ! Test the input parameters. INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KL.LT.0 ) THEN INFO = -3 ELSE IF( KU.LT.0 ) THEN INFO = -4 ELSE IF( LDAB.LT.KL+KV+1 ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGBTRF', -INFO ) RETURN END IF ! Quick return if possible IF( M.EQ.0 .OR. N.EQ.0 ) RETURN ! Determine the block size for this environment NB = ILAENV( 1, 'CGBTRF', ' ', M, N, KL, KU ) ! The block size must not exceed the limit set by the size of the ! local arrays WORK13 and WORK31. NB = MIN( NB, NBMAX ) IF( NB.LE.1 .OR. NB.GT.KL ) THEN ! Use unblocked code CALL CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) ELSE ! Use blocked code ! Zero the superdiagonal elements of the work array WORK13 DO 20 J = 1, NB DO 10 I = 1, J - 1 WORK13( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ! Zero the subdiagonal elements of the work array WORK31 DO 40 J = 1, NB DO 30 I = J + 1, NB WORK31( I, J ) = ZERO 30 CONTINUE 40 CONTINUE ! Gaussian elimination with partial pivoting ! Set fill-in elements in columns KU+2 to KV to zero DO 60 J = KU + 2, MIN( KV, N ) DO 50 I = KV - J + 2, KL AB( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ! JU is the index of the last column affected by the current ! stage of the factorization JU = 1 DO 180 J = 1, MIN( M, N ), NB JB = MIN( NB, MIN( M, N )-J+1 ) ! The active part of the matrix is partitioned ! A11 A12 A13 ! A21 A22 A23 ! A31 A32 A33 ! Here A11, A21 and A31 denote the current block of JB columns ! which is about to be factorized. The number of rows in the ! partitioning are JB, I2, I3 respectively, and the numbers ! of columns are JB, J2, J3. The superdiagonal elements of A13 ! and the subdiagonal elements of A31 lie outside the band. I2 = MIN( KL-JB, M-J-JB+1 ) I3 = MIN( JB, M-J-KL+1 ) ! J2 and J3 are computed after JU has been updated. ! Factorize the current block of JB columns DO 80 JJ = J, J + JB - 1 ! Set fill-in elements in column JJ+KV to zero IF( JJ+KV.LE.N ) THEN DO 70 I = 1, KL AB( I, JJ+KV ) = ZERO 70 CONTINUE END IF ! Find pivot and test for singularity. KM is the number of ! subdiagonal elements in the current column. KM = MIN( KL, M-JJ ) JP = ICAMAX( KM+1, AB( KV+1, JJ ), 1 ) IPIV( JJ ) = JP + JJ - J IF( AB( KV+JP, JJ ).NE.ZERO ) THEN JU = MAX( JU, MIN( JJ+KU+JP-1, N ) ) IF( JP.NE.1 ) THEN ! Apply interchange to columns J to J+JB-1 ! IF( JP+JJ-1.LT.J+KL ) THEN CALL CSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1, & AB( KV+JP+JJ-J, J ), LDAB-1 ) ELSE ! The interchange affects columns J to JJ-1 of A31 ! which are stored in the work array WORK31 CALL CSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, & WORK31( JP+JJ-J-KL, 1 ), LDWORK ) CALL CSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1, & AB( KV+JP, JJ ), LDAB-1 ) END IF END IF ! Compute multipliers CALL CSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), 1) ! Update trailing submatrix within the band and within ! the current block. JM is the index of the last column ! which needs to be updated. JM = MIN( JU, J+JB-1 ) IF( JM.GT.JJ ) & CALL CGERU( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1, & AB( KV, JJ+1 ), LDAB-1, & AB( KV+1, JJ+1 ), LDAB-1 ) ELSE ! If pivot is zero, set INFO to the index of the pivot ! unless a zero pivot has already been found. IF( INFO.EQ.0 ) INFO = JJ END IF ! Copy current column of A31 into the work array WORK31 ! NW = MIN( JJ-J+1, I3 ) IF( NW.GT.0 ) & CALL CCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1, & WORK31( 1, JJ-J+1 ), 1 ) 80 CONTINUE IF( J+JB.LE.N ) THEN ! Apply the row interchanges to the other blocks. J2 = MIN( JU-J+1, KV ) - JB J3 = MAX( 0, JU-J-KV+1 ) ! Use CLASWP to apply the row interchanges to A12, A22, and A32. CALL CLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB, & IPIV( J ), 1 ) ! Adjust the pivot indices. DO 90 I = J, J + JB - 1 IPIV( I ) = IPIV( I ) + J - 1 90 CONTINUE ! Apply the row interchanges to A13, A23, and A33 ! columnwise. K2 = J - 1 + JB + J2 DO 110 I = 1, J3 JJ = K2 + I DO 100 II = J + I - 1, J + JB - 1 IP = IPIV( II ) IF( IP.NE.II ) THEN TEMP = AB( KV+1+II-JJ, JJ ) AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ ) AB( KV+1+IP-JJ, JJ ) = TEMP END IF 100 CONTINUE 110 CONTINUE ! Update the relevant part of the trailing submatrix IF( J2.GT.0 ) THEN ! Update A12 CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Unit', & JB, J2, ONE, AB( KV+1, J ), LDAB-1, & AB( KV+1-JB, J+JB ), LDAB-1 ) IF( I2.GT.0 ) THEN ! Update A22 CALL CGEMM( 'No transpose', 'No transpose', I2, J2, & JB, -ONE, AB( KV+1+JB, J ), LDAB-1, & AB( KV+1-JB, J+JB ), LDAB-1, ONE, & AB( KV+1, J+JB ), LDAB-1 ) END IF IF( I3.GT.0 ) THEN ! Update A32 CALL CGEMM( 'No transpose', 'No transpose', I3, J2, & JB, -ONE, WORK31, LDWORK, & AB( KV+1-JB, J+JB ), LDAB-1, ONE, & AB( KV+KL+1-JB, J+JB ), LDAB-1 ) END IF END IF IF( J3.GT.0 ) THEN ! Copy the lower triangle of A13 into the work array WORK13 DO 130 JJ = 1, J3 DO 120 II = JJ, JB WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 ) 120 CONTINUE 130 CONTINUE ! Update A13 in the work array CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Unit', & JB, J3, ONE, AB( KV+1, J ), LDAB-1, & WORK13, LDWORK ) IF( I2.GT.0 ) THEN ! Update A23 CALL CGEMM( 'No transpose', 'No transpose', I2, J3, & JB, -ONE, AB( KV+1+JB, J ), LDAB-1, & WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), & LDAB-1 ) END IF IF( I3.GT.0 ) THEN ! Update A33 CALL CGEMM( 'No transpose', 'No transpose', I3, J3, & JB, -ONE, WORK31, LDWORK, WORK13, & LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) END IF ! Copy the lower triangle of A13 back into place DO 150 JJ = 1, J3 DO 140 II = JJ, JB AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ ) 140 CONTINUE 150 CONTINUE END IF ELSE ! Adjust the pivot indices. DO 160 I = J, J + JB - 1 IPIV( I ) = IPIV( I ) + J - 1 160 CONTINUE END IF ! Partially undo the interchanges in the current block to ! restore the upper triangular form of A31 and copy the upper ! triangle of A31 back into place DO 170 JJ = J + JB - 1, J, -1 JP = IPIV( JJ ) - JJ + 1 IF( JP.NE.1 ) THEN ! Apply interchange to columns J to JJ-1 IF( JP+JJ-1.LT.J+KL ) THEN ! The interchange does not affect A31 CALL CSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, & AB( KV+JP+JJ-J, J ), LDAB-1 ) ELSE ! The interchange does affect A31 CALL CSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, & WORK31( JP+JJ-J-KL, 1 ), LDWORK ) END IF END IF ! Copy the current column of A31 back into place NW = MIN( I3, JJ-J+1 ) IF( NW.GT.0 ) & CALL CCOPY( NW, WORK31( 1, JJ-J+1 ), 1, & AB( KV+KL+1-JJ+J, JJ ), 1 ) 170 CONTINUE 180 CONTINUE END IF RETURN END SUBROUTINE CGBTRF SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) USE LAPACK2 USE LAPACK1 ! -- LAPACK routine (version 3.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! September 30, 1994 ! .. Scalar Arguments .. INTEGER INFO, KL, KU, LDAB, M, N ! .. ! .. Array Arguments .. INTEGER IPIV( * ) COMPLEX AB( LDAB, * ) ! Purpose ! ======= ! CGBTF2 computes an LU factorization of a complex m-by-n band matrix ! A using partial pivoting with row interchanges. ! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! Arguments ! ========= ! M (input) INTEGER ! The number of rows of the matrix A. M >= 0. ! N (input) INTEGER ! The number of columns of the matrix A. N >= 0. ! KL (input) INTEGER ! The number of subdiagonals within the band of A. KL >= 0. ! KU (input) INTEGER ! The number of superdiagonals within the band of A. KU >= 0. ! AB (input/output) COMPLEX array, dimension (LDAB,N) ! On entry, the matrix A in band storage, in rows KL+1 to ! 2*KL+KU+1; rows 1 to KL of the array need not be set. ! The j-th column of A is stored in the j-th column of the ! array AB as follows: ! AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) ! On exit, details of the factorization: U is stored as an ! upper triangular band matrix with KL+KU superdiagonals in ! rows 1 to KL+KU+1, and the multipliers used during the ! factorization are stored in rows KL+KU+2 to 2*KL+KU+1. ! See below for further details. ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= 2*KL+KU+1. ! IPIV (output) INTEGER array, dimension (min(M,N)) ! The pivot indices; for 1 <= i <= min(M,N), row i of the ! matrix was interchanged with row IPIV(i). ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! > 0: if INFO = +i, U(i,i) is exactly zero. The factorization ! has been completed, but the factor U is exactly ! singular, and division by zero will occur if it is used ! to solve a system of equations. ! Further Details ! =============== ! The band storage scheme is illustrated by the following example, when ! M = N = 6, KL = 2, KU = 1: ! On entry: On exit: ! * * * + + + * * * u14 u25 u36 ! * * + + + + * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * ! a31 a42 a53 a64 * * m31 m42 m53 m64 * * ! Array elements marked * are not used by the routine; elements marked ! + need not be set on entry, but are required by the routine to store ! elements of U, because of fill-in resulting from the row ! interchanges. ! ====================================================================== ! .. Parameters .. COMPLEX ONE, ZERO PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), & ZERO = ( 0.0E+0, 0.0E+0 ) ) ! .. Local Scalars .. INTEGER I, J, JP, JU, KM, KV ! .. External Functions .. ! INTEGER ICAMAX ! EXTERNAL ICAMAX ! .. External Subroutines .. ! EXTERNAL CGERU, CSCAL, CSWAP, XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. Executable Statements .. ! KV is the number of superdiagonals in the factor U, allowing for ! fill-in. KV = KU + KL ! Test the input parameters. INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KL.LT.0 ) THEN INFO = -3 ELSE IF( KU.LT.0 ) THEN INFO = -4 ELSE IF( LDAB.LT.KL+KV+1 ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGBTF2', -INFO ) RETURN END IF ! Quick return if possible IF( M.EQ.0 .OR. N.EQ.0 ) RETURN ! Gaussian elimination with partial pivoting ! Set fill-in elements in columns KU+2 to KV to zero. DO 20 J = KU + 2, MIN( KV, N ) DO 10 I = KV - J + 2, KL AB( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ! JU is the index of the last column affected by the current stage ! of the factorization. JU = 1 DO 40 J = 1, MIN( M, N ) ! Set fill-in elements in column J+KV to zero. IF( J+KV.LE.N ) THEN DO 30 I = 1, KL AB( I, J+KV ) = ZERO 30 CONTINUE END IF ! Find pivot and test for singularity. KM is the number of ! subdiagonal elements in the current column. KM = MIN( KL, M-J ) JP = ICAMAX( KM+1, AB( KV+1, J ), 1 ) IPIV( J ) = JP + J - 1 IF( AB( KV+JP, J ).NE.ZERO ) THEN JU = MAX( JU, MIN( J+KU+JP-1, N ) ) ! Apply interchange to columns J to JU. IF( JP.NE.1 ) & CALL CSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1, & AB( KV+1, J ), LDAB-1 ) IF( KM.GT.0 ) THEN ! Compute multipliers. CALL CSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 ) ! Update trailing submatrix within the band. IF( JU.GT.J ) & CALL CGERU( KM, JU-J, -ONE, AB( KV+2, J ), 1, & AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ), & LDAB-1 ) END IF ELSE ! If pivot is zero, set INFO to the index of the pivot ! unless a zero pivot has already been found. IF( INFO.EQ.0 ) INFO = J END IF 40 CONTINUE RETURN END SUBROUTINE CGBTF2 SUBROUTINE CGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO ) USE LAPACK2 USE LAPACK1 ! ! -- LAPACK routine (version 3.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! September 30, 1994 ! .. Scalar Arguments .. CHARACTER TRANS INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS ! .. Array Arguments .. INTEGER IPIV( * ) COMPLEX AB( LDAB, * ), B( LDB, * ) ! Purpose ! ======= ! CGBTRS solves a system of linear equations ! A * X = B, A**T * X = B, or A**H * X = B ! with a general band matrix A using the LU factorization computed ! by CGBTRF. ! Arguments ! ========= ! TRANS (input) CHARACTER*1 ! Specifies the form of the system of equations. ! = 'N': A * X = B (No transpose) ! = 'T': A**T * X = B (Transpose) ! = 'C': A**H * X = B (Conjugate transpose) ! N (input) INTEGER ! The order of the matrix A. N >= 0. ! KL (input) INTEGER ! The number of subdiagonals within the band of A. KL >= 0. ! KU (input) INTEGER ! The number of superdiagonals within the band of A. KU >= 0. ! NRHS (input) INTEGER ! The number of right hand sides, i.e., the number of columns ! of the matrix B. NRHS >= 0. ! AB (input) COMPLEX array, dimension (LDAB,N) ! Details of the LU factorization of the band matrix A, as ! computed by CGBTRF. U is stored as an upper triangular band ! matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and ! the multipliers used during the factorization are stored in ! rows KL+KU+2 to 2*KL+KU+1. ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= 2*KL+KU+1. ! IPIV (input) INTEGER array, dimension (N) ! The pivot indices; for 1 <= i <= N, row i of the matrix was ! interchanged with row IPIV(i). ! B (input/output) COMPLEX array, dimension (LDB,NRHS) ! On entry, the right hand side matrix B. ! On exit, the solution matrix X. ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(1,N). ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! ! ====================================================================== ! .. Parameters .. COMPLEX ONE PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) ! .. ! .. Local Scalars .. LOGICAL LNOTI, NOTRAN INTEGER I, J, KD, L, LM ! .. ! .. External Functions .. ! LOGICAL LSAME ! EXTERNAL LSAME ! .. ! .. External Subroutines .. ! EXTERNAL CGEMV, CGERU, CLACGV, CSWAP, CTBSV, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. Executable Statements .. ! Test the input parameters. INFO = 0 NOTRAN = LSAME( TRANS, 'N' ) IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. & LSAME( TRANS, 'C' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KL.LT.0 ) THEN INFO = -3 ELSE IF( KU.LT.0 ) THEN INFO = -4 ELSE IF( NRHS.LT.0 ) THEN INFO = -5 ELSE IF( LDAB.LT.( 2*KL+KU+1 ) ) THEN INFO = -7 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGBTRS', -INFO ) RETURN END IF ! Quick return if possible IF( N.EQ.0 .OR. NRHS.EQ.0 ) RETURN KD = KU + KL + 1 LNOTI = KL.GT.0 IF( NOTRAN ) THEN ! Solve A*X = B. ! Solve L*X = B, overwriting B with X. ! L is represented as a product of permutations and unit lower ! triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1), ! where each transformation L(i) is a rank-one modification of ! the identity matrix. IF( LNOTI ) THEN DO 10 J = 1, N - 1 LM = MIN( KL, N-J ) L = IPIV( J ) IF( L.NE.J ) & CALL CSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB ) CALL CGERU( LM, NRHS, -ONE, AB( KD+1, J ), 1, B( J, 1 ), & LDB, B( J+1, 1 ), LDB ) 10 CONTINUE END IF DO 20 I = 1, NRHS ! Solve U*X = B, overwriting B with X. CALL CTBSV( 'Upper', 'No transpose', 'Non-unit', N, KL+KU, & AB, LDAB, B( 1, I ), 1 ) 20 CONTINUE ELSE IF( LSAME( TRANS, 'T' ) ) THEN ! Solve A**T * X = B. DO 30 I = 1, NRHS ! Solve U**T * X = B, overwriting B with X. CALL CTBSV( 'Upper', 'Transpose', 'Non-unit', N, KL+KU, AB, & LDAB, B( 1, I ), 1 ) 30 CONTINUE ! Solve L**T * X = B, overwriting B with X. IF( LNOTI ) THEN DO 40 J = N - 1, 1, -1 LM = MIN( KL, N-J ) CALL CGEMV( 'Transpose', LM, NRHS, -ONE, B( J+1, 1 ), & LDB, AB( KD+1, J ), 1, ONE, B( J, 1 ), LDB ) L = IPIV( J ) IF( L.NE.J ) & CALL CSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB ) 40 CONTINUE END IF ELSE ! Solve A**H * X = B. DO 50 I = 1, NRHS ! Solve U**H * X = B, overwriting B with X. CALL CTBSV( 'Upper', 'Conjugate transpose', 'Non-unit', N, & KL+KU, AB, LDAB, B( 1, I ), 1 ) 50 CONTINUE ! Solve L**H * X = B, overwriting B with X. IF( LNOTI ) THEN DO 60 J = N - 1, 1, -1 LM = MIN( KL, N-J ) CALL CLACGV( NRHS, B( J, 1 ), LDB ) CALL CGEMV( 'Conjugate transpose', LM, NRHS, -ONE, & B( J+1, 1 ), LDB, AB( KD+1, J ), 1, ONE, & B( J, 1 ), LDB ) CALL CLACGV( NRHS, B( J, 1 ), LDB ) L = IPIV( J ) IF( L.NE.J ) & CALL CSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB ) 60 CONTINUE END IF END IF RETURN END SUBROUTINE CGBTRS SUBROUTINE CTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) USE LAPACK2 USE LAPACK1 ! .. Scalar Arguments .. INTEGER INCX, K, LDA, N CHARACTER*1 DIAG, TRANS, UPLO ! .. Array Arguments .. COMPLEX A( LDA, * ), X( * ) ! Purpose ! ======= ! CTBSV solves one of the systems of equations ! A*x = b, or A'*x = b, or conjg( A' )*x = b, ! where b and x are n element vectors and A is an n by n unit, or ! non-unit, upper or lower triangular band matrix, with ( k + 1 ) ! diagonals. ! No test for singularity or near-singularity is included in this ! routine. Such tests must be performed before calling this routine. ! Parameters ! ========== ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! Unchanged on exit. ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the equations to be solved as ! follows: ! TRANS = 'N' or 'n' A*x = b. ! TRANS = 'T' or 't' A'*x = b. ! TRANS = 'C' or 'c' conjg( A' )*x = b. ! Unchanged on exit. ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit ! triangular as follows: ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! Unchanged on exit. ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! K - INTEGER. ! On entry with UPLO = 'U' or 'u', K specifies the number of ! super-diagonals of the matrix A. ! On entry with UPLO = 'L' or 'l', K specifies the number of ! sub-diagonals of the matrix A. ! K must satisfy 0 .le. K. ! Unchanged on exit. ! A - COMPLEX array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) ! by n part of the array A must contain the upper triangular ! band part of the matrix of coefficients, supplied column by ! column, with the leading diagonal of the matrix in row ! ( k + 1 ) of the array, the first super-diagonal starting at ! position 2 in row k, and so on. The top left k by k triangle ! of the array A is not referenced. ! The following program segment will transfer an upper ! triangular band matrix from conventional full matrix storage ! to band storage: ! DO 20, J = 1, N ! M = K + 1 - J ! DO 10, I = MAX( 1, J - K ), J ! A( M + I, J ) = matrix( I, J ) ! 10 CONTINUE ! 20 CONTINUE ! ! Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) ! by n part of the array A must contain the lower triangular ! band part of the matrix of coefficients, supplied column by ! column, with the leading diagonal of the matrix in row 1 of ! the array, the first sub-diagonal starting at position 1 in ! row 2, and so on. The bottom right k by k triangle of the ! array A is not referenced. ! The following program segment will transfer a lower ! triangular band matrix from conventional full matrix storage ! to band storage: ! DO 20, J = 1, N ! M = 1 - J ! DO 10, I = J, MIN( N, J + K ) ! A( M + I, J ) = matrix( I, J ) ! 10 CONTINUE ! 20 CONTINUE ! Note that when DIAG = 'U' or 'u' the elements of the array A ! corresponding to the diagonal elements of the matrix are not ! referenced, but are assumed to be unity. ! Unchanged on exit. ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! ( k + 1 ). ! Unchanged on exit. ! X - COMPLEX array of dimension at least ! ( 1 + ( n - 1 )*abs( INCX ) ). ! Before entry, the incremented array X must contain the n ! element right-hand side vector b. On exit, X is overwritten ! with the solution vector x. ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ====================================================================== ! .. Parameters .. COMPLEX ZERO PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) ! .. Local Scalars .. COMPLEX TEMP INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L LOGICAL NOCONJ, NOUNIT ! .. External Functions .. ! LOGICAL LSAME ! EXTERNAL LSAME ! .. External Subroutines .. ! EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC CONJG, MAX, MIN ! .. Executable Statements .. ! Test the input parameters. INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. & .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. & .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( K.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 7 ELSE IF( INCX.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'CTBSV ', INFO ) RETURN END IF ! Quick return if possible. IF( N.EQ.0 ) RETURN NOCONJ = LSAME( TRANS, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) ! Set up the start point in X if the increment is not unity. This ! will be ( N - 1 )*INCX too small for descending loops. IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF ! Start the operations. In this version the elements of A are ! accessed by sequentially with one pass through A. IF( LSAME( TRANS, 'N' ) )THEN ! Form x := inv( A )*x. IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN L = KPLUS1 - J IF( NOUNIT ) & X( J ) = X( J )/A( KPLUS1, J ) TEMP = X( J ) DO 10, I = J - 1, MAX( 1, J - K ), -1 X( I ) = X( I ) - TEMP*A( L + I, J ) 10 CONTINUE END IF 20 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 40, J = N, 1, -1 KX = KX - INCX IF( X( JX ).NE.ZERO )THEN IX = KX L = KPLUS1 - J IF( NOUNIT ) & X( JX ) = X( JX )/A( KPLUS1, J ) TEMP = X( JX ) DO 30, I = J - 1, MAX( 1, J - K ), -1 X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX - INCX 30 CONTINUE END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN L = 1 - J IF( NOUNIT ) & X( J ) = X( J )/A( 1, J ) TEMP = X( J ) DO 50, I = J + 1, MIN( N, J + K ) X( I ) = X( I ) - TEMP*A( L + I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N KX = KX + INCX IF( X( JX ).NE.ZERO )THEN IX = KX L = 1 - J IF( NOUNIT ) & X( JX ) = X( JX )/A( 1, J ) TEMP = X( JX ) DO 70, I = J + 1, MIN( N, J + K ) X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX + INCX 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE ! Form x := inv( A' )*x or x := inv( conjg( A') )*x. IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = X( J ) L = KPLUS1 - J IF( NOCONJ )THEN DO 90, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) & TEMP = TEMP/A( KPLUS1, J ) ELSE DO 100, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - CONJG( A( L + I, J ) )*X( I ) 100 CONTINUE IF( NOUNIT ) & TEMP = TEMP/CONJG( A( KPLUS1, J ) ) END IF X( J ) = TEMP 110 CONTINUE ELSE JX = KX DO 140, J = 1, N TEMP = X( JX ) IX = KX L = KPLUS1 - J IF( NOCONJ )THEN DO 120, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX + INCX 120 CONTINUE IF( NOUNIT ) & TEMP = TEMP/A( KPLUS1, J ) ELSE DO 130, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - CONJG( A( L + I, J ) )*X( IX ) IX = IX + INCX 130 CONTINUE IF( NOUNIT ) & TEMP = TEMP/CONJG( A( KPLUS1, J ) ) END IF X( JX ) = TEMP JX = JX + INCX IF( J.GT.K ) KX = KX + INCX 140 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 170, J = N, 1, -1 TEMP = X( J ) L = 1 - J IF( NOCONJ )THEN DO 150, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( I ) 150 CONTINUE IF( NOUNIT ) & TEMP = TEMP/A( 1, J ) ELSE DO 160, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - CONJG( A( L + I, J ) )*X( I ) 160 CONTINUE IF( NOUNIT ) & TEMP = TEMP/CONJG( A( 1, J ) ) END IF X( J ) = TEMP 170 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 200, J = N, 1, -1 TEMP = X( JX ) IX = KX L = 1 - J IF( NOCONJ )THEN DO 180, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX - INCX 180 CONTINUE IF( NOUNIT ) & TEMP = TEMP/A( 1, J ) ELSE DO 190, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - CONJG( A( L + I, J ) )*X( IX ) IX = IX - INCX 190 CONTINUE IF( NOUNIT ) & TEMP = TEMP/CONJG( A( 1, J ) ) END IF X( JX ) = TEMP JX = JX - INCX IF( ( N - J ).GE.K ) KX = KX - INCX 200 CONTINUE END IF END IF END IF RETURN END SUBROUTINE CTBSV SUBROUTINE CGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) USE LAPACK2 ! .. Scalar Arguments .. COMPLEX ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS ! .. Array Arguments .. COMPLEX A( LDA, * ), X( * ), Y( * ) ! Purpose ! ======= ! CGEMV performs one of the matrix-vector operations ! y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or ! y := alpha*conjg( A' )*x + beta*y, ! where alpha and beta are scalars, x and y are vectors and A is an ! m by n matrix. ! Parameters ! ========== ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! TRANS = 'N' or 'n' y := alpha*A*x + beta*y. ! TRANS = 'T' or 't' y := alpha*A'*x + beta*y. ! TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y. ! Unchanged on exit. ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix A. ! M must be at least zero. ! Unchanged on exit. ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ALPHA - COMPLEX . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! A - COMPLEX array of DIMENSION ( LDA, n ). ! Before entry, the leading m by n part of the array A must ! contain the matrix of coefficients. ! Unchanged on exit. ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, m ). ! Unchanged on exit. ! X - COMPLEX array of DIMENSION at least ! ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. ! Before entry, the incremented array X must contain the ! vector x. ! Unchanged on exit. ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! BETA - COMPLEX . ! On entry, BETA specifies the scalar beta. When BETA is ! supplied as zero then Y need not be set on input. ! Unchanged on exit. ! Y - COMPLEX array of DIMENSION at least ! ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. ! Before entry with BETA non-zero, the incremented array Y ! must contain the vector y. On exit, Y is overwritten by the ! updated vector y. ! INCY - INTEGER. ! On entry, INCY specifies the increment for the elements of ! Y. INCY must not be zero. ! Unchanged on exit. ! Level 2 Blas routine. ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ====================================================================== ! .. Parameters .. COMPLEX ONE PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) COMPLEX ZERO PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) ! .. Local Scalars .. COMPLEX TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY LOGICAL NOCONJ ! .. External Functions .. ! LOGICAL LSAME ! EXTERNAL LSAME ! .. External Subroutines .. ! EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC CONJG, MAX ! .. Executable Statements .. ! Test the input parameters. ! 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( 'CGEMV ', INFO ) RETURN END IF ! Quick return if possible. IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. & ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) & RETURN NOCONJ = LSAME( TRANS, 'T' ) ! Set LENX and LENY, the lengths of the vectors x and y, and set ! up the start points in X and Y. 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 ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through A. ! First form y := beta*y. 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 ! Form y := alpha*A*x + y. 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 ! Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y. JY = KY IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = ZERO IF( NOCONJ )THEN DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE ELSE DO 100, I = 1, M TEMP = TEMP + CONJG( A( I, J ) )*X( I ) 100 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 110 CONTINUE ELSE DO 140, J = 1, N TEMP = ZERO IX = KX IF( NOCONJ )THEN DO 120, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 120 CONTINUE ELSE DO 130, I = 1, M TEMP = TEMP + CONJG( A( I, J ) )*X( IX ) IX = IX + INCX 130 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 140 CONTINUE END IF END IF RETURN END SUBROUTINE CGEMV SUBROUTINE CLACGV( N, X, INCX ) ! -- LAPACK auxiliary routine (version 3.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! .. Scalar Arguments .. INTEGER INCX, N ! .. Array Arguments .. COMPLEX X( * ) ! Purpose ! ======= ! CLACGV conjugates a complex vector of length N. ! Arguments ! ========= ! N (input) INTEGER ! The length of the vector X. N >= 0. ! X (input/output) COMPLEX array, dimension ! (1+(N-1)*abs(INCX)) ! On entry, the vector of length N to be conjugated. ! On exit, X is overwritten with conjg(X). ! INCX (input) INTEGER ! The spacing between successive elements of X. ! ===================================================================== ! .. Local Scalars .. INTEGER I, IOFF ! .. ! .. Intrinsic Functions .. INTRINSIC CONJG ! .. Executable Statements .. IF( INCX.EQ.1 ) THEN DO 10 I = 1, N X( I ) = CONJG( X( I ) ) 10 CONTINUE ELSE IOFF = 1 IF( INCX.LT.0 ) & IOFF = 1 - ( N-1 )*INCX DO 20 I = 1, N X( IOFF ) = CONJG( X( IOFF ) ) IOFF = IOFF + INCX 20 CONTINUE END IF RETURN END SUBROUTINE CLACGV subroutine ccopy(n,cx,incx,cy,incy) ! copies a vector, x, to a vector, y. ! jack dongarra, linpack, 3/11/78. ! modified 12/3/93, array(1) declarations changed to array(*) ! =============================================================== complex cx(*),cy(*) integer i,incx,incy,ix,iy,n ! if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 ! code for unequal increments or equal increments ! not equal to 1 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 cy(iy) = cx(ix) ix = ix + incx iy = iy + incy 10 continue return ! code for both increments equal to 1 20 do 30 i = 1,n cy(i) = cx(i) 30 continue return end subroutine ccopy END MODULE LAPACK3 !end of file Lapack3.f90