Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/lapack_extra.f
93a1e9712d04feature/refactoring

Last Change
Raw File

lapack_extra.f
View Options

!>
!> @file lapack_extra.f
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
DOUBLE PRECISION FUNCTION DOPGB( SUBNAM, M, N, KL, KU, IPIV )
*
* -- LAPACK timing routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
CHARACTER*6 SUBNAM
INTEGER KL, KU, M, N
* ..
* .. Array Arguments ..
INTEGER IPIV( * )
* ..
*
* Purpose
* =======
*
* DOPGB counts operations for the LU factorization of a band matrix
* xGBTRF.
*
* Arguments
* =========
*
* SUBNAM (input) CHARACTER*6
* The name of the subroutine.
*
* M (input) INTEGER
* The number of rows of the coefficient matrix. M >= 0.
*
* N (input) INTEGER
* The number of columns of the coefficient matrix. N >= 0.
*
* KL (input) INTEGER
* The number of subdiagonals of the matrix. KL >= 0.
*
* KU (input) INTEGER
* The number of superdiagonals of the matrix. KU >= 0.
*
* IPIV (input) INTEGER array, dimension (min(M,N))
* The vector of pivot indices from DGBTRF or ZGBTRF.
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL CORZ, SORD
CHARACTER C1
CHARACTER*2 C2
CHARACTER*3 C3
INTEGER I, J, JP, JU, KM
DOUBLE PRECISION ADDFAC, ADDS, MULFAC, MULTS
* ..
* .. External Functions ..
LOGICAL LSAME, LSAMEN
EXTERNAL LSAME, LSAMEN
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
DOPGB = 0
MULTS = 0
ADDS = 0
C1 = SUBNAM( 1: 1 )
C2 = SUBNAM( 2: 3 )
C3 = SUBNAM( 4: 6 )
SORD = LSAME( C1, 'S' ) .OR. LSAME( C1, 'D' )
CORZ = LSAME( C1, 'C' ) .OR. LSAME( C1, 'Z' )
IF( .NOT.( SORD .OR. CORZ ) )
$ RETURN
IF( LSAME( C1, 'S' ) .OR. LSAME( C1, 'D' ) ) THEN
ADDFAC = 1
MULFAC = 1
ELSE
ADDFAC = 2
MULFAC = 6
END IF
*
* --------------------------
* GB: General Band matrices
* --------------------------
*
IF( LSAMEN( 2, C2, 'GB' ) ) THEN
*
* xGBTRF: M, N, KL, KU => M, N, KL, KU
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
JU = 1
DO 10 J = 1, MIN( M, N )
KM = MIN( KL, M-J )
JP = IPIV( J )
JU = MAX( JU, MIN( JP+KU, N ) )
IF( KM.GT.0 ) THEN
MULTS = MULTS + KM*( 1+JU-J )
ADDS = ADDS + KM*( JU-J )
END IF
10 CONTINUE
END IF
*
* ---------------------------------
* GT: General Tridiagonal matrices
* ---------------------------------
*
ELSE IF( LSAMEN( 2, C2, 'GT' ) ) THEN
*
* xGTTRF: N => M
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
MULTS = 2*( M-1 )
ADDS = M - 1
DO 20 I = 1, M - 2
IF( IPIV( I ).NE.I )
$ MULTS = MULTS + 1
20 CONTINUE
*
* xGTTRS: N, NRHS => M, N
*
ELSE IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = 4*N*( M-1 )
ADDS = 3*N*( M-1 )
*
* xGTSV: N, NRHS => M, N
*
ELSE IF( LSAMEN( 3, C3, 'SV ' ) ) THEN
MULTS = ( 4*N+2 )*( M-1 )
ADDS = ( 3*N+1 )*( M-1 )
DO 30 I = 1, M - 2
IF( IPIV( I ).NE.I )
$ MULTS = MULTS + 1
30 CONTINUE
END IF
END IF
*
DOPGB = MULFAC*MULTS + ADDFAC*ADDS
RETURN
*
* End of DOPGB
*
END
DOUBLE PRECISION FUNCTION DOPLA( SUBNAM, M, N, KL, KU, NB )
*
* -- LAPACK timing routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
CHARACTER*6 SUBNAM
INTEGER KL, KU, M, N, NB
* ..
*
* Purpose
* =======
*
* DOPLA computes an approximation of the number of floating point
* operations used by the subroutine SUBNAM with the given values
* of the parameters M, N, KL, KU, and NB.
*
* This version counts operations for the LAPACK subroutines.
*
* Arguments
* =========
*
* SUBNAM (input) CHARACTER*6
* The name of the subroutine.
*
* M (input) INTEGER
* The number of rows of the coefficient matrix. M >= 0.
*
* N (input) INTEGER
* The number of columns of the coefficient matrix.
* For solve routine when the matrix is square,
* N is the number of right hand sides. N >= 0.
*
* KL (input) INTEGER
* The lower band width of the coefficient matrix.
* If needed, 0 <= KL <= M-1.
* For xGEQRS, KL is the number of right hand sides.
*
* KU (input) INTEGER
* The upper band width of the coefficient matrix.
* If needed, 0 <= KU <= N-1.
*
* NB (input) INTEGER
* The block size. If needed, NB >= 1.
*
* Notes
* =====
*
* In the comments below, the association is given between arguments
* in the requested subroutine and local arguments. For example,
*
* xGETRS: N, NRHS => M, N
*
* means that arguments N and NRHS in DGETRS are passed to arguments
* M and N in this procedure.
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL CORZ, SORD
CHARACTER C1
CHARACTER*2 C2
CHARACTER*3 C3
INTEGER I
DOUBLE PRECISION ADDFAC, ADDS, EK, EM, EMN, EN, MULFAC, MULTS,
$ WL, WU
* ..
* .. External Functions ..
LOGICAL LSAME, LSAMEN
EXTERNAL LSAME, LSAMEN
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
* --------------------------------------------------------
* Initialize DOPLA to 0 and do a quick return if possible.
* --------------------------------------------------------
*
DOPLA = 0
MULTS = 0
ADDS = 0
C1 = SUBNAM( 1: 1 )
C2 = SUBNAM( 2: 3 )
C3 = SUBNAM( 4: 6 )
SORD = LSAME( C1, 'S' ) .OR. LSAME( C1, 'D' )
CORZ = LSAME( C1, 'C' ) .OR. LSAME( C1, 'Z' )
IF( M.LE.0 .OR. .NOT.( SORD .OR. CORZ ) )
$ RETURN
*
* ---------------------------------------------------------
* If the coefficient matrix is real, count each add as 1
* operation and each multiply as 1 operation.
* If the coefficient matrix is complex, count each add as 2
* operations and each multiply as 6 operations.
* ---------------------------------------------------------
*
IF( LSAME( C1, 'S' ) .OR. LSAME( C1, 'D' ) ) THEN
ADDFAC = 1
MULFAC = 1
ELSE
ADDFAC = 2
MULFAC = 6
END IF
EM = M
EN = N
EK = KL
*
* ---------------------------------
* GE: GEneral rectangular matrices
* ---------------------------------
*
IF( LSAMEN( 2, C2, 'GE' ) ) THEN
*
* xGETRF: M, N => M, N
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
EMN = MIN( M, N )
ADDS = EMN*( EM*EN-( EM+EN )*( EMN+1.D0 ) / 2.D0+
$ ( EMN+1.D0 )*( 2.D0*EMN+1.D0 ) / 6.D0 )
MULTS = ADDS + EMN*( EM-( EMN+1.D0 ) / 2.D0 )
*
* xGETRS: N, NRHS => M, N
*
ELSE IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = EN*EM*EM
ADDS = EN*( EM*( EM-1.D0 ) )
*
* xGETRI: N => M
*
ELSE IF( LSAMEN( 3, C3, 'TRI' ) ) THEN
MULTS = EM*( 5.D0 / 6.D0+EM*( 1.D0 / 2.D0+EM*( 2.D0 /
$ 3.D0 ) ) )
ADDS = EM*( 5.D0 / 6.D0+EM*( -3.D0 / 2.D0+EM*( 2.D0 /
$ 3.D0 ) ) )
*
* xGEQRF or xGEQLF: M, N => M, N
*
ELSE IF( LSAMEN( 3, C3, 'QRF' ) .OR.
$ LSAMEN( 3, C3, 'QR2' ) .OR.
$ LSAMEN( 3, C3, 'QLF' ) .OR. LSAMEN( 3, C3, 'QL2' ) )
$ THEN
IF( M.GE.N ) THEN
MULTS = EN*( ( ( 23.D0 / 6.D0 )+EM+EN / 2.D0 )+EN*
$ ( EM-EN / 3.D0 ) )
ADDS = EN*( ( 5.D0 / 6.D0 )+EN*
$ ( 1.D0 / 2.D0+( EM-EN / 3.D0 ) ) )
ELSE
MULTS = EM*( ( ( 23.D0 / 6.D0 )+2.D0*EN-EM / 2.D0 )+EM*
$ ( EN-EM / 3.D0 ) )
ADDS = EM*( ( 5.D0 / 6.D0 )+EN-EM / 2.D0+EM*
$ ( EN-EM / 3.D0 ) )
END IF
*
* xGERQF or xGELQF: M, N => M, N
*
ELSE IF( LSAMEN( 3, C3, 'RQF' ) .OR.
$ LSAMEN( 3, C3, 'RQ2' ) .OR.
$ LSAMEN( 3, C3, 'LQF' ) .OR. LSAMEN( 3, C3, 'LQ2' ) )
$ THEN
IF( M.GE.N ) THEN
MULTS = EN*( ( ( 29.D0 / 6.D0 )+EM+EN / 2.D0 )+EN*
$ ( EM-EN / 3.D0 ) )
ADDS = EN*( ( 5.D0 / 6.D0 )+EM+EN*
$ ( -1.D0 / 2.D0+( EM-EN / 3.D0 ) ) )
ELSE
MULTS = EM*( ( ( 29.D0 / 6.D0 )+2.D0*EN-EM / 2.D0 )+EM*
$ ( EN-EM / 3.D0 ) )
ADDS = EM*( ( 5.D0 / 6.D0 )+EM / 2.D0+EM*
$ ( EN-EM / 3.D0 ) )
END IF
*
* xGEQPF: M, N => M, N
*
ELSE IF( LSAMEN( 3, C3, 'QPF' ) ) THEN
EMN = MIN( M, N )
MULTS = 2*EN*EN + EMN*( 3*EM+5*EN+2*EM*EN-( EMN+1 )*
$ ( 4+EN+EM-( 2*EMN+1 ) / 3 ) )
ADDS = EN*EN + EMN*( 2*EM+EN+2*EM*EN-( EMN+1 )*
$ ( 2+EN+EM-( 2*EMN+1 ) / 3 ) )
*
* xGEQRS or xGERQS: M, N, NRHS => M, N, KL
*
ELSE IF( LSAMEN( 3, C3, 'QRS' ) .OR. LSAMEN( 3, C3, 'RQS' ) )
$ THEN
MULTS = EK*( EN*( 2.D0-EK )+EM*
$ ( 2.D0*EN+( EM+1.D0 ) / 2.D0 ) )
ADDS = EK*( EN*( 1.D0-EK )+EM*
$ ( 2.D0*EN+( EM-1.D0 ) / 2.D0 ) )
*
* xGELQS or xGEQLS: M, N, NRHS => M, N, KL
*
ELSE IF( LSAMEN( 3, C3, 'LQS' ) .OR. LSAMEN( 3, C3, 'QLS' ) )
$ THEN
MULTS = EK*( EM*( 2.D0-EK )+EN*
$ ( 2.D0*EM+( EN+1.D0 ) / 2.D0 ) )
ADDS = EK*( EM*( 1.D0-EK )+EN*
$ ( 2.D0*EM+( EN-1.D0 ) / 2.D0 ) )
*
* xGEBRD: M, N => M, N
*
ELSE IF( LSAMEN( 3, C3, 'BRD' ) ) THEN
IF( M.GE.N ) THEN
MULTS = EN*( 20.D0 / 3.D0+EN*
$ ( 2.D0+( 2.D0*EM-( 2.D0 / 3.D0 )*EN ) ) )
ADDS = EN*( 5.D0 / 3.D0+( EN-EM )+EN*
$ ( 2.D0*EM-( 2.D0 / 3.D0 )*EN ) )
ELSE
MULTS = EM*( 20.D0 / 3.D0+EM*
$ ( 2.D0+( 2.D0*EN-( 2.D0 / 3.D0 )*EM ) ) )
ADDS = EM*( 5.D0 / 3.D0+( EM-EN )+EM*
$ ( 2.D0*EN-( 2.D0 / 3.D0 )*EM ) )
END IF
*
* xGEHRD: N => M
*
ELSE IF( LSAMEN( 3, C3, 'HRD' ) ) THEN
IF( M.EQ.1 ) THEN
MULTS = 0.D0
ADDS = 0.D0
ELSE
MULTS = -13.D0 + EM*( -7.D0 / 6.D0+EM*
$ ( 0.5D0+EM*( 5.D0 / 3.D0 ) ) )
ADDS = -8.D0 + EM*( -2.D0 / 3.D0+EM*
$ ( -1.D0+EM*( 5.D0 / 3.D0 ) ) )
END IF
*
END IF
*
* ----------------------------
* GB: General Banded matrices
* ----------------------------
* Note: The operation count is overestimated because
* it is assumed that the factor U fills in to the maximum
* extent, i.e., that its bandwidth goes from KU to KL + KU.
*
ELSE IF( LSAMEN( 2, C2, 'GB' ) ) THEN
*
* xGBTRF: M, N, KL, KU => M, N, KL, KU
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
DO 10 I = MIN( M, N ), 1, -1
WL = MAX( 0, MIN( KL, M-I ) )
WU = MAX( 0, MIN( KL+KU, N-I ) )
MULTS = MULTS + WL*( 1.D0+WU )
ADDS = ADDS + WL*WU
10 CONTINUE
*
* xGBTRS: N, NRHS, KL, KU => M, N, KL, KU
*
ELSE IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
WL = MAX( 0, MIN( KL, M-1 ) )
WU = MAX( 0, MIN( KL+KU, M-1 ) )
MULTS = EN*( EM*( WL+1.D0+WU )-0.5D0*
$ ( WL*( WL+1.D0 )+WU*( WU+1.D0 ) ) )
ADDS = EN*( EM*( WL+WU )-0.5D0*
$ ( WL*( WL+1.D0 )+WU*( WU+1.D0 ) ) )
*
END IF
*
* --------------------------------------
* PO: POsitive definite matrices
* PP: Positive definite Packed matrices
* --------------------------------------
*
ELSE IF( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'PP' ) ) THEN
*
* xPOTRF: N => M
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
MULTS = EM*( 1.D0 / 3.D0+EM*( 1.D0 / 2.D0+EM*( 1.D0 /
$ 6.D0 ) ) )
ADDS = ( 1.D0 / 6.D0 )*EM*( -1.D0+EM*EM )
*
* xPOTRS: N, NRHS => M, N
*
ELSE IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = EN*( EM*( EM+1.D0 ) )
ADDS = EN*( EM*( EM-1.D0 ) )
*
* xPOTRI: N => M
*
ELSE IF( LSAMEN( 3, C3, 'TRI' ) ) THEN
MULTS = EM*( 2.D0 / 3.D0+EM*( 1.D0+EM*( 1.D0 / 3.D0 ) ) )
ADDS = EM*( 1.D0 / 6.D0+EM*( -1.D0 / 2.D0+EM*( 1.D0 /
$ 3.D0 ) ) )
*
END IF
*
* ------------------------------------
* PB: Positive definite Band matrices
* ------------------------------------
*
ELSE IF( LSAMEN( 2, C2, 'PB' ) ) THEN
*
* xPBTRF: N, K => M, KL
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
MULTS = EK*( -2.D0 / 3.D0+EK*( -1.D0+EK*( -1.D0 / 3.D0 ) ) )
$ + EM*( 1.D0+EK*( 3.D0 / 2.D0+EK*( 1.D0 / 2.D0 ) ) )
ADDS = EK*( -1.D0 / 6.D0+EK*( -1.D0 / 2.D0+EK*( -1.D0 /
$ 3.D0 ) ) ) + EM*( EK / 2.D0*( 1.D0+EK ) )
*
* xPBTRS: N, NRHS, K => M, N, KL
*
ELSE IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = EN*( ( 2*EM-EK )*( EK+1.D0 ) )
ADDS = EN*( EK*( 2*EM-( EK+1.D0 ) ) )
*
END IF
*
* ----------------------------------
* PT: Positive definite Tridiagonal
* ----------------------------------
*
ELSE IF( LSAMEN( 2, C2, 'PT' ) ) THEN
*
* xPTTRF: N => M
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
MULTS = 2*( EM-1 )
ADDS = EM - 1
*
* xPTTRS: N, NRHS => M, N
*
ELSE IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = EN*( 3*EM-2 )
ADDS = EN*( 2*( EM-1 ) )
*
* xPTSV: N, NRHS => M, N
*
ELSE IF( LSAMEN( 3, C3, 'SV ' ) ) THEN
MULTS = 2*( EM-1 ) + EN*( 3*EM-2 )
ADDS = EM - 1 + EN*( 2*( EM-1 ) )
END IF
*
* --------------------------------------------------------
* SY: SYmmetric indefinite matrices
* SP: Symmetric indefinite Packed matrices
* HE: HErmitian indefinite matrices (complex only)
* HP: Hermitian indefinite Packed matrices (complex only)
* --------------------------------------------------------
*
ELSE IF( LSAMEN( 2, C2, 'SY' ) .OR. LSAMEN( 2, C2, 'SP' ) .OR.
$ LSAMEN( 3, SUBNAM, 'CHE' ) .OR.
$ LSAMEN( 3, SUBNAM, 'ZHE' ) .OR.
$ LSAMEN( 3, SUBNAM, 'CHP' ) .OR.
$ LSAMEN( 3, SUBNAM, 'ZHP' ) ) THEN
*
* xSYTRF: N => M
*
IF( LSAMEN( 3, C3, 'TRF' ) ) THEN
MULTS = EM*( 10.D0 / 3.D0+EM*
$ ( 1.D0 / 2.D0+EM*( 1.D0 / 6.D0 ) ) )
ADDS = EM / 6.D0*( -1.D0+EM*EM )
*
* xSYTRS: N, NRHS => M, N
*
ELSE IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = EN*EM*EM
ADDS = EN*( EM*( EM-1.D0 ) )
*
* xSYTRI: N => M
*
ELSE IF( LSAMEN( 3, C3, 'TRI' ) ) THEN
MULTS = EM*( 2.D0 / 3.D0+EM*EM*( 1.D0 / 3.D0 ) )
ADDS = EM*( -1.D0 / 3.D0+EM*EM*( 1.D0 / 3.D0 ) )
*
* xSYTRD, xSYTD2: N => M
*
ELSE IF( LSAMEN( 3, C3, 'TRD' ) .OR. LSAMEN( 3, C3, 'TD2' ) )
$ THEN
IF( M.EQ.1 ) THEN
MULTS = 0.D0
ADDS = 0.D0
ELSE
MULTS = -15.D0 + EM*( -1.D0 / 6.D0+EM*
$ ( 5.D0 / 2.D0+EM*( 2.D0 / 3.D0 ) ) )
ADDS = -4.D0 + EM*( -8.D0 / 3.D0+EM*
$ ( 1.D0+EM*( 2.D0 / 3.D0 ) ) )
END IF
END IF
*
* -------------------
* Triangular matrices
* -------------------
*
ELSE IF( LSAMEN( 2, C2, 'TR' ) .OR. LSAMEN( 2, C2, 'TP' ) ) THEN
*
* xTRTRS: N, NRHS => M, N
*
IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = EN*EM*( EM+1.D0 ) / 2.D0
ADDS = EN*EM*( EM-1.D0 ) / 2.D0
*
* xTRTRI: N => M
*
ELSE IF( LSAMEN( 3, C3, 'TRI' ) ) THEN
MULTS = EM*( 1.D0 / 3.D0+EM*( 1.D0 / 2.D0+EM*( 1.D0 /
$ 6.D0 ) ) )
ADDS = EM*( 1.D0 / 3.D0+EM*( -1.D0 / 2.D0+EM*( 1.D0 /
$ 6.D0 ) ) )
*
END IF
*
ELSE IF( LSAMEN( 2, C2, 'TB' ) ) THEN
*
* xTBTRS: N, NRHS, K => M, N, KL
*
IF( LSAMEN( 3, C3, 'TRS' ) ) THEN
MULTS = EN*( EM*( EM+1.D0 ) / 2.D0-( EM-EK-1.D0 )*
$ ( EM-EK ) / 2.D0 )
ADDS = EN*( EM*( EM-1.D0 ) / 2.D0-( EM-EK-1.D0 )*( EM-EK ) /
$ 2.D0 )
END IF
*
* --------------------
* Trapezoidal matrices
* --------------------
*
ELSE IF( LSAMEN( 2, C2, 'TZ' ) ) THEN
*
* xTZRQF: M, N => M, N
*
IF( LSAMEN( 3, C3, 'RQF' ) ) THEN
EMN = MIN( M, N )
MULTS = 3*EM*( EN-EM+1 ) + ( 2*EN-2*EM+3 )*
$ ( EM*EM-EMN*( EMN+1 ) / 2 )
ADDS = ( EN-EM+1 )*( EM+2*EM*EM-EMN*( EMN+1 ) )
END IF
*
* -------------------
* Orthogonal matrices
* -------------------
*
ELSE IF( ( SORD .AND. LSAMEN( 2, C2, 'OR' ) ) .OR.
$ ( CORZ .AND. LSAMEN( 2, C2, 'UN' ) ) ) THEN
*
* -MQR, -MLQ, -MQL, or -MRQ: M, N, K, SIDE => M, N, KL, KU
* where KU<= 0 indicates SIDE = 'L'
* and KU> 0 indicates SIDE = 'R'
*
IF( LSAMEN( 3, C3, 'MQR' ) .OR. LSAMEN( 3, C3, 'MLQ' ) .OR.
$ LSAMEN( 3, C3, 'MQL' ) .OR. LSAMEN( 3, C3, 'MRQ' ) ) THEN
IF( KU.LE.0 ) THEN
MULTS = EK*EN*( 2.D0*EM+2.D0-EK )
ADDS = EK*EN*( 2.D0*EM+1.D0-EK )
ELSE
MULTS = EK*( EM*( 2.D0*EN-EK )+
$ ( EM+EN+( 1.D0-EK ) / 2.D0 ) )
ADDS = EK*EM*( 2.D0*EN+1.D0-EK )
END IF
*
* -GQR or -GQL: M, N, K => M, N, KL
*
ELSE IF( LSAMEN( 3, C3, 'GQR' ) .OR. LSAMEN( 3, C3, 'GQL' ) )
$ THEN
MULTS = EK*( -5.D0 / 3.D0+( 2.D0*EN-EK )+
$ ( 2.D0*EM*EN+EK*( ( 2.D0 / 3.D0 )*EK-EM-EN ) ) )
ADDS = EK*( 1.D0 / 3.D0+( EN-EM )+
$ ( 2.D0*EM*EN+EK*( ( 2.D0 / 3.D0 )*EK-EM-EN ) ) )
*
* -GLQ or -GRQ: M, N, K => M, N, KL
*
ELSE IF( LSAMEN( 3, C3, 'GLQ' ) .OR. LSAMEN( 3, C3, 'GRQ' ) )
$ THEN
MULTS = EK*( -2.D0 / 3.D0+( EM+EN-EK )+
$ ( 2.D0*EM*EN+EK*( ( 2.D0 / 3.D0 )*EK-EM-EN ) ) )
ADDS = EK*( 1.D0 / 3.D0+( EM-EN )+
$ ( 2.D0*EM*EN+EK*( ( 2.D0 / 3.D0 )*EK-EM-EN ) ) )
*
END IF
*
END IF
*
DOPLA = MULFAC*MULTS + ADDFAC*ADDS
*
RETURN
*
* End of DOPLA
*
END
LOGICAL FUNCTION LSAMEN( N, CA, CB )
*
* -- LAPACK auxiliary routine (version 2.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
CHARACTER*( * ) CA, CB
INTEGER N
* ..
*
* Purpose
* =======
*
* LSAMEN tests if the first N letters of CA are the same as the
* first N letters of CB, regardless of case.
* LSAMEN returns .TRUE. if CA and CB are equivalent except for case
* and .FALSE. otherwise. LSAMEN also returns .FALSE. if LEN( CA )
* or LEN( CB ) is less than N.
*
* Arguments
* =========
*
* N (input) INTEGER
* The number of characters in CA and CB to be compared.
*
* CA (input) CHARACTER*(*)
* CB (input) CHARACTER*(*)
* CA and CB specify two character strings of length at least N.
* Only the first N characters of each string will be accessed.
*
* =====================================================================
*
* .. Local Scalars ..
INTEGER I
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. Intrinsic Functions ..
INTRINSIC LEN
* ..
* .. Executable Statements ..
*
LSAMEN = .FALSE.
IF( LEN( CA ).LT.N .OR. LEN( CB ).LT.N )
$ GO TO 20
*
* Do for each character in the two strings.
*
DO 10 I = 1, N
*
* Test if the characters are equal using LSAME.
*
IF( .NOT.LSAME( CA( I: I ), CB( I: I ) ) )
$ GO TO 20
*
10 CONTINUE
LSAMEN = .TRUE.
*
20 CONTINUE
RETURN
*
* End of LSAMEN
*
END
c4science · Help