From b6c0693ab456c51930eebc3634660506e1a3e45b Mon Sep 17 00:00:00 2001 From: Vincent Vanlaer Date: Thu, 8 Aug 2024 00:38:04 -0400 Subject: [PATCH 1/3] mtx: remove lapack quad No code in MESA makes use of this --- mtx/lapack_quad/iqamax.f | 53 -- mtx/lapack_quad/qgemm.f | 313 ------- mtx/lapack_quad/qgemv.f | 261 ------ mtx/lapack_quad/qger.f | 159 ---- mtx/lapack_quad/qgetf2.f | 147 --- mtx/lapack_quad/qgetrf.f | 159 ---- mtx/lapack_quad/qgetrs.f | 149 --- mtx/lapack_quad/qgttrf.f | 168 ---- mtx/lapack_quad/qgttrs.f | 140 --- mtx/lapack_quad/qgtts2.f | 196 ---- mtx/lapack_quad/qlamch.f | 857 ------------------ mtx/lapack_quad/qlaswp.f | 119 --- mtx/lapack_quad/qscal.f | 57 -- mtx/lapack_quad/qswap.f | 70 -- mtx/lapack_quad/qtrsm.f | 373 -------- mtx/lapack_quad/qtrsv.f | 338 ------- mtx/make/makefile_base | 36 +- mtx/private/bcyclic.f90 | 2 +- mtx/private/mtx_support.f90 | 2 +- .../{my_lapack95.F90 => my_lapack95.f90} | 17 +- mtx/public/mtx_lapack95.dek | 6 +- mtx/test/make/makefile_base | 2 +- mtx/test/src/test_mtx.f90 | 2 - mtx/test/src/test_mtx_support.f90 | 40 - mtx/test/src/test_square_quad.f90 | 175 ---- mtx/test/test_output | 15 - 26 files changed, 10 insertions(+), 3846 deletions(-) delete mode 100644 mtx/lapack_quad/iqamax.f delete mode 100644 mtx/lapack_quad/qgemm.f delete mode 100644 mtx/lapack_quad/qgemv.f delete mode 100644 mtx/lapack_quad/qger.f delete mode 100644 mtx/lapack_quad/qgetf2.f delete mode 100644 mtx/lapack_quad/qgetrf.f delete mode 100644 mtx/lapack_quad/qgetrs.f delete mode 100644 mtx/lapack_quad/qgttrf.f delete mode 100644 mtx/lapack_quad/qgttrs.f delete mode 100644 mtx/lapack_quad/qgtts2.f delete mode 100644 mtx/lapack_quad/qlamch.f delete mode 100644 mtx/lapack_quad/qlaswp.f delete mode 100644 mtx/lapack_quad/qscal.f delete mode 100644 mtx/lapack_quad/qswap.f delete mode 100644 mtx/lapack_quad/qtrsm.f delete mode 100644 mtx/lapack_quad/qtrsv.f rename mtx/private/{my_lapack95.F90 => my_lapack95.f90} (99%) delete mode 100644 mtx/test/src/test_square_quad.f90 diff --git a/mtx/lapack_quad/iqamax.f b/mtx/lapack_quad/iqamax.f deleted file mode 100644 index 06c6a8658..000000000 --- a/mtx/lapack_quad/iqamax.f +++ /dev/null @@ -1,53 +0,0 @@ - INTEGER FUNCTION IQAMAX(N,DX,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N -* .. -* .. Array Arguments .. - REAL(16) DX(*) -* .. -* -* Purpose -* ======= -* -* finds the index of element having max. absolute value. -* jack dongarra, linpack, 3/11/78. -* modified 3/93 to return if incx .le. 0. -* modified 12/3/93, array(1) declarations changed to array(*) -* -* -* .. Local Scalars .. - REAL(16) QMAX - INTEGER I,IX -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS -* .. - IQAMAX = 0 - IF (N.LT.1 .OR. INCX.LE.0) RETURN - IQAMAX = 1 - IF (N.EQ.1) RETURN - IF (INCX.EQ.1) GO TO 20 -* -* code for increment not equal to 1 -* - IX = 1 - QMAX = ABS(DX(1)) - IX = IX + INCX - DO 10 I = 2,N - IF (ABS(DX(IX)).LE.QMAX) GO TO 5 - IQAMAX = I - QMAX = ABS(DX(IX)) - 5 IX = IX + INCX - 10 CONTINUE - RETURN -* -* code for increment equal to 1 -* - 20 QMAX = ABS(DX(1)) - DO 30 I = 2,N - IF (ABS(DX(I)).LE.QMAX) GO TO 30 - IQAMAX = I - QMAX = ABS(DX(I)) - 30 CONTINUE - RETURN - END diff --git a/mtx/lapack_quad/qgemm.f b/mtx/lapack_quad/qgemm.f deleted file mode 100644 index 18f5e62e9..000000000 --- a/mtx/lapack_quad/qgemm.f +++ /dev/null @@ -1,313 +0,0 @@ - SUBROUTINE QGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) -* .. Scalar Arguments .. - REAL(16) ALPHA,BETA - INTEGER K,LDA,LDB,LDC,M,N - CHARACTER TRANSA,TRANSB -* .. -* .. Array Arguments .. - REAL(16) A(LDA,*),B(LDB,*),C(LDC,*) -* .. -* -* Purpose -* ======= -* -* QGEMM performs one of the matrix-matrix operations -* -* C := alpha*op( A )*op( B ) + beta*C, -* -* where op( X ) is one of -* -* op( X ) = X or op( X ) = X', -* -* alpha and beta are scalars, and A, B and C are matrices, with op( A ) -* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. -* -* Arguments -* ========== -* -* TRANSA - CHARACTER*1. -* On entry, TRANSA specifies the form of op( A ) to be used in -* the matrix multiplication as follows: -* -* TRANSA = 'N' or 'n', op( A ) = A. -* -* TRANSA = 'T' or 't', op( A ) = A'. -* -* TRANSA = 'C' or 'c', op( A ) = A'. -* -* Unchanged on exit. -* -* TRANSB - CHARACTER*1. -* On entry, TRANSB specifies the form of op( B ) to be used in -* the matrix multiplication as follows: -* -* TRANSB = 'N' or 'n', op( B ) = B. -* -* TRANSB = 'T' or 't', op( B ) = B'. -* -* TRANSB = 'C' or 'c', op( B ) = B'. -* -* Unchanged on exit. -* -* M - INTEGER. -* On entry, M specifies the number of rows of the matrix -* op( A ) and of the matrix C. M must be at least zero. -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the number of columns of the matrix -* op( B ) and the number of columns of the matrix C. N must be -* at least zero. -* Unchanged on exit. -* -* K - INTEGER. -* On entry, K specifies the number of columns of the matrix -* op( A ) and the number of rows of the matrix op( B ). K must -* be at least zero. -* Unchanged on exit. -* -* ALPHA - REAL(16). -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* A - REAL(16) array of DIMENSION ( LDA, ka ), where ka is -* k when TRANSA = 'N' or 'n', and is m otherwise. -* Before entry with TRANSA = 'N' or 'n', the leading m by k -* part of the array A must contain the matrix A, otherwise -* the leading k by m part of the array A must contain the -* matrix A. -* Unchanged on exit. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. When TRANSA = 'N' or 'n' then -* LDA must be at least max( 1, m ), otherwise LDA must be at -* least max( 1, k ). -* Unchanged on exit. -* -* B - REAL(16) array of DIMENSION ( LDB, kb ), where kb is -* n when TRANSB = 'N' or 'n', and is k otherwise. -* Before entry with TRANSB = 'N' or 'n', the leading k by n -* part of the array B must contain the matrix B, otherwise -* the leading n by k part of the array B must contain the -* matrix B. -* Unchanged on exit. -* -* LDB - INTEGER. -* On entry, LDB specifies the first dimension of B as declared -* in the calling (sub) program. When TRANSB = 'N' or 'n' then -* LDB must be at least max( 1, k ), otherwise LDB must be at -* least max( 1, n ). -* Unchanged on exit. -* -* BETA - REAL(16). -* On entry, BETA specifies the scalar beta. When BETA is -* supplied as zero then C need not be set on input. -* Unchanged on exit. -* -* C - REAL(16) array of DIMENSION ( LDC, n ). -* Before entry, the leading m by n part of the array C must -* contain the matrix C, except when beta is zero, in which -* case C need not be set on entry. -* On exit, the array C is overwritten by the m by n matrix -* ( alpha*op( A )*op( B ) + beta*C ). -* -* LDC - INTEGER. -* On entry, LDC specifies the first dimension of C as declared -* in the calling (sub) program. LDC must be at least -* max( 1, m ). -* Unchanged on exit. -* -* -* Level 3 Blas routine. -* -* -- Written on 8-February-1989. -* Jack Dongarra, Argonne National Laboratory. -* Iain Duff, AERE Harwell. -* Jeremy Du Croz, Numerical Algorithms Group Ltd. -* Sven Hammarling, Numerical Algorithms Group Ltd. -* -* -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* .. Local Scalars .. - REAL(16) TEMP - INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB - LOGICAL NOTA,NOTB -* .. -* .. Parameters .. - REAL(16) ONE,ZERO - PARAMETER (ONE=1.0_16,ZERO=0.0_16) -* .. -* -* Set NOTA and NOTB as true if A and B respectively are not -* transposed and set NROWA, NCOLA and NROWB as the number of rows -* and columns of A and the number of rows of B respectively. -* - NOTA = LSAME(TRANSA,'N') - NOTB = LSAME(TRANSB,'N') - IF (NOTA) THEN - NROWA = M - NCOLA = K - ELSE - NROWA = K - NCOLA = M - END IF - IF (NOTB) THEN - NROWB = K - ELSE - NROWB = N - END IF -* -* Test the input parameters. -* - INFO = 0 - IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND. - + (.NOT.LSAME(TRANSA,'T'))) THEN - INFO = 1 - ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND. - + (.NOT.LSAME(TRANSB,'T'))) THEN - INFO = 2 - ELSE IF (M.LT.0) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (K.LT.0) THEN - INFO = 5 - ELSE IF (LDA.LT.MAX(1,NROWA)) THEN - INFO = 8 - ELSE IF (LDB.LT.MAX(1,NROWB)) THEN - INFO = 10 - ELSE IF (LDC.LT.MAX(1,M)) THEN - INFO = 13 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('QGEMM ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((M.EQ.0) .OR. (N.EQ.0) .OR. - + (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN -* -* And if alpha.eq.zero. -* - IF (ALPHA.EQ.ZERO) THEN - IF (BETA.EQ.ZERO) THEN - DO 20 J = 1,N - DO 10 I = 1,M - C(I,J) = ZERO - 10 CONTINUE - 20 CONTINUE - ELSE - DO 40 J = 1,N - DO 30 I = 1,M - C(I,J) = BETA*C(I,J) - 30 CONTINUE - 40 CONTINUE - END IF - RETURN - END IF -* -* Start the operations. -* - IF (NOTB) THEN - IF (NOTA) THEN -* -* Form C := alpha*A*B + beta*C. -* - DO 90 J = 1,N - IF (BETA.EQ.ZERO) THEN - DO 50 I = 1,M - C(I,J) = ZERO - 50 CONTINUE - ELSE IF (BETA.NE.ONE) THEN - DO 60 I = 1,M - C(I,J) = BETA*C(I,J) - 60 CONTINUE - END IF - DO 80 L = 1,K - IF (B(L,J).NE.ZERO) THEN - TEMP = ALPHA*B(L,J) - DO 70 I = 1,M - C(I,J) = C(I,J) + TEMP*A(I,L) - 70 CONTINUE - END IF - 80 CONTINUE - 90 CONTINUE - ELSE -* -* Form C := alpha*A'*B + beta*C -* - DO 120 J = 1,N - DO 110 I = 1,M - TEMP = ZERO - DO 100 L = 1,K - TEMP = TEMP + A(L,I)*B(L,J) - 100 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 110 CONTINUE - 120 CONTINUE - END IF - ELSE - IF (NOTA) THEN -* -* Form C := alpha*A*B' + beta*C -* - DO 170 J = 1,N - IF (BETA.EQ.ZERO) THEN - DO 130 I = 1,M - C(I,J) = ZERO - 130 CONTINUE - ELSE IF (BETA.NE.ONE) THEN - DO 140 I = 1,M - C(I,J) = BETA*C(I,J) - 140 CONTINUE - END IF - DO 160 L = 1,K - IF (B(J,L).NE.ZERO) THEN - TEMP = ALPHA*B(J,L) - DO 150 I = 1,M - C(I,J) = C(I,J) + TEMP*A(I,L) - 150 CONTINUE - END IF - 160 CONTINUE - 170 CONTINUE - ELSE -* -* Form C := alpha*A'*B' + beta*C -* - DO 200 J = 1,N - DO 190 I = 1,M - TEMP = ZERO - DO 180 L = 1,K - TEMP = TEMP + A(L,I)*B(J,L) - 180 CONTINUE - IF (BETA.EQ.ZERO) THEN - C(I,J) = ALPHA*TEMP - ELSE - C(I,J) = ALPHA*TEMP + BETA*C(I,J) - END IF - 190 CONTINUE - 200 CONTINUE - END IF - END IF -* - RETURN -* -* End of QGEMM . -* - END diff --git a/mtx/lapack_quad/qgemv.f b/mtx/lapack_quad/qgemv.f deleted file mode 100644 index 7f655ce5b..000000000 --- a/mtx/lapack_quad/qgemv.f +++ /dev/null @@ -1,261 +0,0 @@ - SUBROUTINE QGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) -* .. Scalar Arguments .. - real(16) ALPHA,BETA - INTEGER INCX,INCY,LDA,M,N - CHARACTER TRANS -* .. -* .. Array Arguments .. - real(16) A(LDA,*),X(*),Y(*) -* .. -* -* Purpose -* ======= -* -* QGEMV performs one of the matrix-vector operations -* -* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, -* -* where alpha and beta are scalars, x and y are vectors and A is an -* m by n matrix. -* -* Arguments -* ========== -* -* 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*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 - real(16). -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* A - real(16) 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 - real(16) 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 - real(16). -* 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 - real(16) 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 .. - real(16) ONE,ZERO - PARAMETER (ONE=1.0_16,ZERO=0.0_16) -* .. -* .. Local Scalars .. - real(16) TEMP - INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* -* 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('QGEMV ',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 -* -* 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. -* - JY = KY - IF (INCX.EQ.1) THEN - DO 100 J = 1,N - TEMP = ZERO - DO 90 I = 1,M - TEMP = TEMP + A(I,J)*X(I) - 90 CONTINUE - Y(JY) = Y(JY) + ALPHA*TEMP - JY = JY + INCY - 100 CONTINUE - ELSE - DO 120 J = 1,N - TEMP = ZERO - IX = KX - DO 110 I = 1,M - TEMP = TEMP + A(I,J)*X(IX) - IX = IX + INCX - 110 CONTINUE - Y(JY) = Y(JY) + ALPHA*TEMP - JY = JY + INCY - 120 CONTINUE - END IF - END IF -* - RETURN -* -* End of QGEMV . -* - END diff --git a/mtx/lapack_quad/qger.f b/mtx/lapack_quad/qger.f deleted file mode 100644 index a9b246fe9..000000000 --- a/mtx/lapack_quad/qger.f +++ /dev/null @@ -1,159 +0,0 @@ - SUBROUTINE QGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) -* .. Scalar Arguments .. - REAL(16) ALPHA - INTEGER INCX,INCY,LDA,M,N -* .. -* .. Array Arguments .. - REAL(16) A(LDA,*),X(*),Y(*) -* .. -* -* Purpose -* ======= -* -* QGER performs the rank 1 operation -* -* A := alpha*x*y' + A, -* -* where alpha is a scalar, x is an m element vector, y is an n element -* vector and A is an m by n matrix. -* -* Arguments -* ========== -* -* 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 - REAL(16). -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - REAL(16) array of dimension at least -* ( 1 + ( m - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the m -* element 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. -* -* Y - REAL(16) array of dimension at least -* ( 1 + ( n - 1 )*abs( INCY ) ). -* Before entry, the incremented array Y must contain the n -* element vector y. -* Unchanged on exit. -* -* INCY - INTEGER. -* On entry, INCY specifies the increment for the elements of -* Y. INCY must not be zero. -* Unchanged on exit. -* -* A - REAL(16) array of DIMENSION ( LDA, n ). -* Before entry, the leading m by n part of the array A must -* contain the matrix of coefficients. On exit, A is -* overwritten by the updated matrix. -* -* 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. -* -* -* 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 .. - REAL(16) ZERO - PARAMETER (ZERO=0.0_16) -* .. -* .. Local Scalars .. - REAL(16) TEMP - INTEGER I,INFO,IX,J,JY,KX -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (M.LT.0) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - ELSE IF (INCY.EQ.0) THEN - INFO = 7 - ELSE IF (LDA.LT.MAX(1,M)) THEN - INFO = 9 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('QGER ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((M.EQ.0) .OR. (N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Start the operations. In this version the elements of A are -* accessed sequentially with one pass through A. -* - IF (INCY.GT.0) THEN - JY = 1 - ELSE - JY = 1 - (N-1)*INCY - END IF - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (Y(JY).NE.ZERO) THEN - TEMP = ALPHA*Y(JY) - DO 10 I = 1,M - A(I,J) = A(I,J) + X(I)*TEMP - 10 CONTINUE - END IF - JY = JY + INCY - 20 CONTINUE - ELSE - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (M-1)*INCX - END IF - DO 40 J = 1,N - IF (Y(JY).NE.ZERO) THEN - TEMP = ALPHA*Y(JY) - IX = KX - DO 30 I = 1,M - A(I,J) = A(I,J) + X(IX)*TEMP - IX = IX + INCX - 30 CONTINUE - END IF - JY = JY + INCY - 40 CONTINUE - END IF -* - RETURN -* -* End of QGER . -* - END diff --git a/mtx/lapack_quad/qgetf2.f b/mtx/lapack_quad/qgetf2.f deleted file mode 100644 index c124a3dba..000000000 --- a/mtx/lapack_quad/qgetf2.f +++ /dev/null @@ -1,147 +0,0 @@ - SUBROUTINE QGETF2( M, N, A, LDA, IPIV, INFO ) -* -* -- LAPACK routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. Scalar Arguments .. - INTEGER INFO, LDA, M, N -* .. -* .. Array Arguments .. - INTEGER IPIV( * ) - REAL(16) A( LDA, * ) -* .. -* -* Purpose -* ======= -* -* QGETF2 computes an LU factorization of a general m-by-n matrix A -* using partial pivoting with row interchanges. -* -* The factorization has the form -* A = P * L * U -* where P is a permutation matrix, L is lower triangular with unit -* diagonal elements (lower trapezoidal if m > n), and U is upper -* triangular (upper trapezoidal if m < n). -* -* This is the right-looking Level 2 BLAS version of the algorithm. -* -* 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. -* -* A (input/output) REAL(16) array, dimension (LDA,N) -* On entry, the m by n matrix to be factored. -* On exit, the factors L and U from the factorization -* A = P*L*U; the unit diagonal elements of L are not stored. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,M). -* -* 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 = -k, the k-th argument had an illegal value -* > 0: if INFO = k, U(k,k) 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. -* -* ===================================================================== -* -* .. Parameters .. - REAL(16) ONE, ZERO - PARAMETER ( ONE = 1.0_16, ZERO = 0.0_16 ) -* .. -* .. Local Scalars .. - REAL(16) SFMIN - INTEGER I, J, JP -* .. -* .. External Functions .. - REAL(16) QLAMCH - INTEGER IQAMAX - EXTERNAL QLAMCH, IQAMAX -* .. -* .. External Subroutines .. - EXTERNAL QGER, QSCAL, QSWAP, XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX, MIN -* .. -* .. Executable Statements .. -* -* Test the input parameters. -* - INFO = 0 - IF( M.LT.0 ) THEN - INFO = -1 - ELSE IF( N.LT.0 ) THEN - INFO = -2 - ELSE IF( LDA.LT.MAX( 1, M ) ) THEN - INFO = -4 - END IF - IF( INFO.NE.0 ) THEN - CALL XERBLA( 'QGETF2', -INFO ) - RETURN - END IF -* -* Quick return if possible -* - IF( M.EQ.0 .OR. N.EQ.0 ) - $ RETURN -* -* Compute machine safe minimum -* - SFMIN = QLAMCH('S') -* - DO 10 J = 1, MIN( M, N ) -* -* Find pivot and test for singularity. -* - JP = J - 1 + IQAMAX( M-J+1, A( J, J ), 1 ) - IPIV( J ) = JP - IF( A( JP, J ).NE.ZERO ) THEN -* -* Apply the interchange to columns 1:N. -* - IF( JP.NE.J ) - $ CALL QSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) -* -* Compute elements J+1:M of J-th column. -* - IF( J.LT.M ) THEN - IF( ABS(A( J, J )) .GE. SFMIN ) THEN - CALL QSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) - ELSE - DO 20 I = 1, M-J - A( J+I, J ) = A( J+I, J ) / A( J, J ) - 20 CONTINUE - END IF - END IF -* - ELSE IF( INFO.EQ.0 ) THEN -* - INFO = J - END IF -* - IF( J.LT.MIN( M, N ) ) THEN -* -* Update trailing submatrix. -* - CALL QGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, - $ A( J+1, J+1 ), LDA ) - END IF - 10 CONTINUE - RETURN -* -* End of QGETF2 -* - END diff --git a/mtx/lapack_quad/qgetrf.f b/mtx/lapack_quad/qgetrf.f deleted file mode 100644 index f16de3b5e..000000000 --- a/mtx/lapack_quad/qgetrf.f +++ /dev/null @@ -1,159 +0,0 @@ - SUBROUTINE QGETRF( M, N, A, LDA, IPIV, INFO ) -* -* -- LAPACK routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. Scalar Arguments .. - INTEGER INFO, LDA, M, N -* .. -* .. Array Arguments .. - INTEGER IPIV( * ) - REAL(16) A( LDA, * ) -* .. -* -* Purpose -* ======= -* -* QGETRF computes an LU factorization of a general M-by-N matrix A -* using partial pivoting with row interchanges. -* -* The factorization has the form -* A = P * L * U -* where P is a permutation matrix, L is lower triangular with unit -* diagonal elements (lower trapezoidal if m > n), and U is upper -* triangular (upper trapezoidal if m < n). -* -* This is the right-looking Level 3 BLAS version of the algorithm. -* -* 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. -* -* A (input/output) REAL(16) array, dimension (LDA,N) -* On entry, the M-by-N matrix to be factored. -* On exit, the factors L and U from the factorization -* A = P*L*U; the unit diagonal elements of L are not stored. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,M). -* -* 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. -* -* ===================================================================== -* -* .. Parameters .. - REAL(16) ONE - PARAMETER ( ONE = 1.0_16 ) -* .. -* .. Local Scalars .. - INTEGER I, IINFO, J, JB, NB -* .. -* .. External Subroutines .. - EXTERNAL QGEMM, QGETF2, QLASWP, QTRSM, XERBLA -* .. -* .. External Functions .. - INTEGER ILAENV - EXTERNAL ILAENV -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX, MIN -* .. -* .. Executable Statements .. -* -* Test the input parameters. -* - INFO = 0 - IF( M.LT.0 ) THEN - INFO = -1 - ELSE IF( N.LT.0 ) THEN - INFO = -2 - ELSE IF( LDA.LT.MAX( 1, M ) ) THEN - INFO = -4 - END IF - IF( INFO.NE.0 ) THEN - CALL XERBLA( 'QGETRF', -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, 'QGETRF', ' ', M, N, -1, -1 ) - IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN -* -* Use unblocked code. -* - CALL QGETF2( M, N, A, LDA, IPIV, INFO ) - ELSE -* -* Use blocked code. -* - DO 20 J = 1, MIN( M, N ), NB - JB = MIN( MIN( M, N )-J+1, NB ) -* -* Factor diagonal and subdiagonal blocks and test for exact -* singularity. -* - CALL QGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO ) -* -* Adjust INFO and the pivot indices. -* - IF( INFO.EQ.0 .AND. IINFO.GT.0 ) - $ INFO = IINFO + J - 1 - DO 10 I = J, MIN( M, J+JB-1 ) - IPIV( I ) = J - 1 + IPIV( I ) - 10 CONTINUE -* -* Apply interchanges to columns 1:J-1. -* - CALL QLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) -* - IF( J+JB.LE.N ) THEN -* -* Apply interchanges to columns J+JB:N. -* - CALL QLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, - $ IPIV, 1 ) -* -* Compute block row of U. -* - CALL QTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, - $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), - $ LDA ) - IF( J+JB.LE.M ) THEN -* -* Update trailing submatrix. -* - CALL QGEMM( 'No transpose', 'No transpose', M-J-JB+1, - $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, - $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), - $ LDA ) - END IF - END IF - 20 CONTINUE - END IF - RETURN -* -* End of QGETRF -* - END diff --git a/mtx/lapack_quad/qgetrs.f b/mtx/lapack_quad/qgetrs.f deleted file mode 100644 index f27645973..000000000 --- a/mtx/lapack_quad/qgetrs.f +++ /dev/null @@ -1,149 +0,0 @@ - SUBROUTINE QGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) -* -* -- LAPACK routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. Scalar Arguments .. - CHARACTER TRANS - INTEGER INFO, LDA, LDB, N, NRHS -* .. -* .. Array Arguments .. - INTEGER IPIV( * ) - REAL(16) A( LDA, * ), B( LDB, * ) -* .. -* -* Purpose -* ======= -* -* QGETRS solves a system of linear equations -* A * X = B or A' * X = B -* with a general N-by-N matrix A using the LU factorization computed -* by DGETRF. -* -* Arguments -* ========= -* -* TRANS (input) CHARACTER*1 -* Specifies the form of the system of equations: -* = 'N': A * X = B (No transpose) -* = 'T': A'* X = B (Transpose) -* = 'C': A'* X = B (Conjugate transpose = Transpose) -* -* N (input) INTEGER -* The order of the matrix A. N >= 0. -* -* NRHS (input) INTEGER -* The number of right hand sides, i.e., the number of columns -* of the matrix B. NRHS >= 0. -* -* A (input) REAL(16) array, dimension (LDA,N) -* The factors L and U from the factorization A = P*L*U -* as computed by DGETRF. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* IPIV (input) INTEGER array, dimension (N) -* The pivot indices from DGETRF; for 1<=i<=N, row i of the -* matrix was interchanged with row IPIV(i). -* -* B (input/output) REAL(16) 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 .. - REAL(16) ONE - PARAMETER ( ONE = 1.0_16 ) -* .. -* .. Local Scalars .. - LOGICAL NOTRAN -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL QLASWP, QTRSM, XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* .. 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( NRHS.LT.0 ) THEN - INFO = -3 - ELSE IF( LDA.LT.MAX( 1, N ) ) THEN - INFO = -5 - ELSE IF( LDB.LT.MAX( 1, N ) ) THEN - INFO = -8 - END IF - IF( INFO.NE.0 ) THEN - CALL XERBLA( 'QGETRS', -INFO ) - RETURN - END IF -* -* Quick return if possible -* - IF( N.EQ.0 .OR. NRHS.EQ.0 ) - $ RETURN -* - IF( NOTRAN ) THEN -* -* Solve A * X = B. -* -* Apply row interchanges to the right hand sides. -* - CALL QLASWP( NRHS, B, LDB, 1, N, IPIV, 1 ) -* -* Solve L*X = B, overwriting B with X. -* - CALL QTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, - $ ONE, A, LDA, B, LDB ) -* -* Solve U*X = B, overwriting B with X. -* - CALL QTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, - $ NRHS, ONE, A, LDA, B, LDB ) - ELSE -* -* Solve A' * X = B. -* -* Solve U'*X = B, overwriting B with X. -* - CALL QTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, - $ ONE, A, LDA, B, LDB ) -* -* Solve L'*X = B, overwriting B with X. -* - CALL QTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE, - $ A, LDA, B, LDB ) -* -* Apply row interchanges to the solution vectors. -* - CALL QLASWP( NRHS, B, LDB, 1, N, IPIV, -1 ) - END IF -* - RETURN -* -* End of QGETRS -* - END diff --git a/mtx/lapack_quad/qgttrf.f b/mtx/lapack_quad/qgttrf.f deleted file mode 100644 index 7480aa6d9..000000000 --- a/mtx/lapack_quad/qgttrf.f +++ /dev/null @@ -1,168 +0,0 @@ - SUBROUTINE QGTTRF( N, DL, D, DU, DU2, IPIV, INFO ) -* -* -- LAPACK routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. Scalar Arguments .. - INTEGER INFO, N -* .. -* .. Array Arguments .. - INTEGER IPIV( * ) - REAL(16) D( * ), DL( * ), DU( * ), DU2( * ) -* .. -* -* Purpose -* ======= -* -* QGTTRF computes an LU factorization of a real tridiagonal matrix A -* using elimination with partial pivoting and row interchanges. -* -* The factorization has the form -* A = L * U -* where L is a product of permutation and unit lower bidiagonal -* matrices and U is upper triangular with nonzeros in only the main -* diagonal and first two superdiagonals. -* -* Arguments -* ========= -* -* N (input) INTEGER -* The order of the matrix A. -* -* DL (input/output) REAL(16) array, dimension (N-1) -* On entry, DL must contain the (n-1) sub-diagonal elements of -* A. -* -* On exit, DL is overwritten by the (n-1) multipliers that -* define the matrix L from the LU factorization of A. -* -* D (input/output) REAL(16) array, dimension (N) -* On entry, D must contain the diagonal elements of A. -* -* On exit, D is overwritten by the n diagonal elements of the -* upper triangular matrix U from the LU factorization of A. -* -* DU (input/output) REAL(16) array, dimension (N-1) -* On entry, DU must contain the (n-1) super-diagonal elements -* of A. -* -* On exit, DU is overwritten by the (n-1) elements of the first -* super-diagonal of U. -* -* DU2 (output) REAL(16) array, dimension (N-2) -* On exit, DU2 is overwritten by the (n-2) elements of the -* second super-diagonal of U. -* -* IPIV (output) INTEGER array, dimension (N) -* The pivot indices; for 1 <= i <= n, row i of the matrix was -* interchanged with row IPIV(i). IPIV(i) will always be either -* i or i+1; IPIV(i) = i indicates a row interchange was not -* required. -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -k, the k-th argument had an illegal value -* > 0: if INFO = k, U(k,k) 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. -* -* ===================================================================== -* -* .. Parameters .. - REAL(16) ZERO - PARAMETER ( ZERO = 0.0_16 ) -* .. -* .. Local Scalars .. - INTEGER I - REAL(16) FACT, TEMP -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Executable Statements .. -* - INFO = 0 - IF( N.LT.0 ) THEN - INFO = -1 - CALL XERBLA( 'QGTTRF', -INFO ) - RETURN - END IF -* -* Quick return if possible -* - IF( N.EQ.0 ) - $ RETURN -* -* Initialize IPIV(i) = i and DU2(I) = 0 -* - DO 10 I = 1, N - IPIV( I ) = I - 10 CONTINUE - DO 20 I = 1, N - 2 - DU2( I ) = ZERO - 20 CONTINUE -* - DO 30 I = 1, N - 2 - IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN -* -* No row interchange required, eliminate DL(I) -* - IF( D( I ).NE.ZERO ) THEN - FACT = DL( I ) / D( I ) - DL( I ) = FACT - D( I+1 ) = D( I+1 ) - FACT*DU( I ) - END IF - ELSE -* -* Interchange rows I and I+1, eliminate DL(I) -* - FACT = D( I ) / DL( I ) - D( I ) = DL( I ) - DL( I ) = FACT - TEMP = DU( I ) - DU( I ) = D( I+1 ) - D( I+1 ) = TEMP - FACT*D( I+1 ) - DU2( I ) = DU( I+1 ) - DU( I+1 ) = -FACT*DU( I+1 ) - IPIV( I ) = I + 1 - END IF - 30 CONTINUE - IF( N.GT.1 ) THEN - I = N - 1 - IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN - IF( D( I ).NE.ZERO ) THEN - FACT = DL( I ) / D( I ) - DL( I ) = FACT - D( I+1 ) = D( I+1 ) - FACT*DU( I ) - END IF - ELSE - FACT = D( I ) / DL( I ) - D( I ) = DL( I ) - DL( I ) = FACT - TEMP = DU( I ) - DU( I ) = D( I+1 ) - D( I+1 ) = TEMP - FACT*D( I+1 ) - IPIV( I ) = I + 1 - END IF - END IF -* -* Check for a zero on the diagonal of U. -* - DO 40 I = 1, N - IF( D( I ).EQ.ZERO ) THEN - INFO = I - GO TO 50 - END IF - 40 CONTINUE - 50 CONTINUE -* - RETURN -* -* End of QGTTRF -* - END diff --git a/mtx/lapack_quad/qgttrs.f b/mtx/lapack_quad/qgttrs.f deleted file mode 100644 index 2b1189763..000000000 --- a/mtx/lapack_quad/qgttrs.f +++ /dev/null @@ -1,140 +0,0 @@ - SUBROUTINE QGTTRS( TRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB, - $ INFO ) -* -* -- LAPACK routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. Scalar Arguments .. - CHARACTER TRANS - INTEGER INFO, LDB, N, NRHS -* .. -* .. Array Arguments .. - INTEGER IPIV( * ) - REAL(16) B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * ) -* .. -* -* Purpose -* ======= -* -* QGTTRS solves one of the systems of equations -* A*X = B or A'*X = B, -* with a tridiagonal matrix A using the LU factorization computed -* by QGTTRF. -* -* Arguments -* ========= -* -* TRANS (input) CHARACTER*1 -* Specifies the form of the system of equations. -* = 'N': A * X = B (No transpose) -* = 'T': A'* X = B (Transpose) -* = 'C': A'* X = B (Conjugate transpose = Transpose) -* -* N (input) INTEGER -* The order of the matrix A. -* -* NRHS (input) INTEGER -* The number of right hand sides, i.e., the number of columns -* of the matrix B. NRHS >= 0. -* -* DL (input) REAL(16) array, dimension (N-1) -* The (n-1) multipliers that define the matrix L from the -* LU factorization of A. -* -* D (input) REAL(16) array, dimension (N) -* The n diagonal elements of the upper triangular matrix U from -* the LU factorization of A. -* -* DU (input) REAL(16) array, dimension (N-1) -* The (n-1) elements of the first super-diagonal of U. -* -* DU2 (input) REAL(16) array, dimension (N-2) -* The (n-2) elements of the second super-diagonal of U. -* -* IPIV (input) INTEGER array, dimension (N) -* The pivot indices; for 1 <= i <= n, row i of the matrix was -* interchanged with row IPIV(i). IPIV(i) will always be either -* i or i+1; IPIV(i) = i indicates a row interchange was not -* required. -* -* B (input/output) REAL(16) array, dimension (LDB,NRHS) -* On entry, the matrix of right hand side vectors B. -* On exit, B is overwritten by the solution vectors 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 -* -* ===================================================================== -* -* .. Local Scalars .. - LOGICAL NOTRAN - INTEGER ITRANS, J, JB, NB -* .. -* .. External Functions .. - INTEGER ILAENV - EXTERNAL ILAENV -* .. -* .. External Subroutines .. - EXTERNAL QGTTS2, XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX, MIN -* .. -* .. Executable Statements .. -* - INFO = 0 - NOTRAN = ( TRANS.EQ.'N' .OR. TRANS.EQ.'n' ) - IF( .NOT.NOTRAN .AND. .NOT.( TRANS.EQ.'T' .OR. TRANS.EQ. - $ 't' ) .AND. .NOT.( TRANS.EQ.'C' .OR. TRANS.EQ.'c' ) ) THEN - INFO = -1 - ELSE IF( N.LT.0 ) THEN - INFO = -2 - ELSE IF( NRHS.LT.0 ) THEN - INFO = -3 - ELSE IF( LDB.LT.MAX( N, 1 ) ) THEN - INFO = -10 - END IF - IF( INFO.NE.0 ) THEN - CALL XERBLA( 'QGTTRS', -INFO ) - RETURN - END IF -* -* Quick return if possible -* - IF( N.EQ.0 .OR. NRHS.EQ.0 ) - $ RETURN -* -* Decode TRANS -* - IF( NOTRAN ) THEN - ITRANS = 0 - ELSE - ITRANS = 1 - END IF -* -* Determine the number of right-hand sides to solve at a time. -* - IF( NRHS.EQ.1 ) THEN - NB = 1 - ELSE - NB = MAX( 1, ILAENV( 1, 'QGTTRS', TRANS, N, NRHS, -1, -1 ) ) - END IF -* - IF( NB.GE.NRHS ) THEN - CALL QGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB ) - ELSE - DO 10 J = 1, NRHS, NB - JB = MIN( NRHS-J+1, NB ) - CALL QGTTS2( ITRANS, N, JB, DL, D, DU, DU2, IPIV, B( 1, J ), - $ LDB ) - 10 CONTINUE - END IF -* -* End of QGTTRS -* - END diff --git a/mtx/lapack_quad/qgtts2.f b/mtx/lapack_quad/qgtts2.f deleted file mode 100644 index 3523521f3..000000000 --- a/mtx/lapack_quad/qgtts2.f +++ /dev/null @@ -1,196 +0,0 @@ - SUBROUTINE QGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB ) -* -* -- LAPACK auxiliary routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. Scalar Arguments .. - INTEGER ITRANS, LDB, N, NRHS -* .. -* .. Array Arguments .. - INTEGER IPIV( * ) - REAL(16) B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * ) -* .. -* -* Purpose -* ======= -* -* QGTTS2 solves one of the systems of equations -* A*X = B or A'*X = B, -* with a tridiagonal matrix A using the LU factorization computed -* by QGTTRF. -* -* Arguments -* ========= -* -* ITRANS (input) INTEGER -* Specifies the form of the system of equations. -* = 0: A * X = B (No transpose) -* = 1: A'* X = B (Transpose) -* = 2: A'* X = B (Conjugate transpose = Transpose) -* -* N (input) INTEGER -* The order of the matrix A. -* -* NRHS (input) INTEGER -* The number of right hand sides, i.e., the number of columns -* of the matrix B. NRHS >= 0. -* -* DL (input) REAL(16) array, dimension (N-1) -* The (n-1) multipliers that define the matrix L from the -* LU factorization of A. -* -* D (input) REAL(16) array, dimension (N) -* The n diagonal elements of the upper triangular matrix U from -* the LU factorization of A. -* -* DU (input) REAL(16) array, dimension (N-1) -* The (n-1) elements of the first super-diagonal of U. -* -* DU2 (input) REAL(16) array, dimension (N-2) -* The (n-2) elements of the second super-diagonal of U. -* -* IPIV (input) INTEGER array, dimension (N) -* The pivot indices; for 1 <= i <= n, row i of the matrix was -* interchanged with row IPIV(i). IPIV(i) will always be either -* i or i+1; IPIV(i) = i indicates a row interchange was not -* required. -* -* B (input/output) REAL(16) array, dimension (LDB,NRHS) -* On entry, the matrix of right hand side vectors B. -* On exit, B is overwritten by the solution vectors X. -* -* LDB (input) INTEGER -* The leading dimension of the array B. LDB >= max(1,N). -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I, IP, J - REAL(16) TEMP -* .. -* .. Executable Statements .. -* -* Quick return if possible -* - IF( N.EQ.0 .OR. NRHS.EQ.0 ) - $ RETURN -* - IF( ITRANS.EQ.0 ) THEN -* -* Solve A*X = B using the LU factorization of A, -* overwriting each right hand side vector with its solution. -* - IF( NRHS.LE.1 ) THEN - J = 1 - 10 CONTINUE -* -* Solve L*x = b. -* - DO 20 I = 1, N - 1 - IP = IPIV( I ) - TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J ) - B( I, J ) = B( IP, J ) - B( I+1, J ) = TEMP - 20 CONTINUE -* -* Solve U*x = b. -* - B( N, J ) = B( N, J ) / D( N ) - IF( N.GT.1 ) - $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / - $ D( N-1 ) - DO 30 I = N - 2, 1, -1 - B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* - $ B( I+2, J ) ) / D( I ) - 30 CONTINUE - IF( J.LT.NRHS ) THEN - J = J + 1 - GO TO 10 - END IF - ELSE - DO 60 J = 1, NRHS -* -* Solve L*x = b. -* - DO 40 I = 1, N - 1 - IF( IPIV( I ).EQ.I ) THEN - B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J ) - ELSE - TEMP = B( I, J ) - B( I, J ) = B( I+1, J ) - B( I+1, J ) = TEMP - DL( I )*B( I, J ) - END IF - 40 CONTINUE -* -* Solve U*x = b. -* - B( N, J ) = B( N, J ) / D( N ) - IF( N.GT.1 ) - $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / - $ D( N-1 ) - DO 50 I = N - 2, 1, -1 - B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* - $ B( I+2, J ) ) / D( I ) - 50 CONTINUE - 60 CONTINUE - END IF - ELSE -* -* Solve A' * X = B. -* - IF( NRHS.LE.1 ) THEN -* -* Solve U'*x = b. -* - J = 1 - 70 CONTINUE - B( 1, J ) = B( 1, J ) / D( 1 ) - IF( N.GT.1 ) - $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) - DO 80 I = 3, N - B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )* - $ B( I-2, J ) ) / D( I ) - 80 CONTINUE -* -* Solve L'*x = b. -* - DO 90 I = N - 1, 1, -1 - IP = IPIV( I ) - TEMP = B( I, J ) - DL( I )*B( I+1, J ) - B( I, J ) = B( IP, J ) - B( IP, J ) = TEMP - 90 CONTINUE - IF( J.LT.NRHS ) THEN - J = J + 1 - GO TO 70 - END IF -* - ELSE - DO 120 J = 1, NRHS -* -* Solve U'*x = b. -* - B( 1, J ) = B( 1, J ) / D( 1 ) - IF( N.GT.1 ) - $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) - DO 100 I = 3, N - B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )- - $ DU2( I-2 )*B( I-2, J ) ) / D( I ) - 100 CONTINUE - DO 110 I = N - 1, 1, -1 - IF( IPIV( I ).EQ.I ) THEN - B( I, J ) = B( I, J ) - DL( I )*B( I+1, J ) - ELSE - TEMP = B( I+1, J ) - B( I+1, J ) = B( I, J ) - DL( I )*TEMP - B( I, J ) = TEMP - END IF - 110 CONTINUE - 120 CONTINUE - END IF - END IF -* -* End of QGTTS2 -* - END diff --git a/mtx/lapack_quad/qlamch.f b/mtx/lapack_quad/qlamch.f deleted file mode 100644 index 671fc42d5..000000000 --- a/mtx/lapack_quad/qlamch.f +++ /dev/null @@ -1,857 +0,0 @@ - REAL(16) FUNCTION QLAMCH( CMACH ) -* -* -- 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 .. - CHARACTER CMACH -* .. -* -* Purpose -* ======= -* -* QLAMCH determines REAL(16) machine parameters. -* -* Arguments -* ========= -* -* CMACH (input) CHARACTER*1 -* Specifies the value to be returned by QLAMCH: -* = 'E' or 'e', QLAMCH := eps -* = 'S' or 's , QLAMCH := sfmin -* = 'B' or 'b', QLAMCH := base -* = 'P' or 'p', QLAMCH := eps*base -* = 'N' or 'n', QLAMCH := t -* = 'R' or 'r', QLAMCH := rnd -* = 'M' or 'm', QLAMCH := emin -* = 'U' or 'u', QLAMCH := rmin -* = 'L' or 'l', QLAMCH := emax -* = 'O' or 'o', QLAMCH := rmax -* -* where -* -* eps = relative machine precision -* sfmin = safe minimum, such that 1/sfmin does not overflow -* base = base of the machine -* prec = eps*base -* t = number of (base) digits in the mantissa -* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise -* emin = minimum exponent before (gradual) underflow -* rmin = underflow threshold - base**(emin-1) -* emax = largest exponent before overflow -* rmax = overflow threshold - (base**emax)*(1-eps) -* -* ===================================================================== -* -* .. Parameters .. - REAL(16) ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) -* .. -* .. Local Scalars .. - LOGICAL FIRST, LRND - INTEGER BETA, IMAX, IMIN, IT - REAL(16) BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, - $ RND, SFMIN, SMALL, T -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL QLAMC2 -* .. -* .. Save statement .. - SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, - $ EMAX, RMAX, PREC -* .. -* .. Data statements .. - DATA FIRST / .TRUE. / -* .. -* .. Executable Statements .. -* - IF( FIRST ) THEN - FIRST = .FALSE. - CALL QLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) - BASE = BETA - T = IT - IF( LRND ) THEN - RND = ONE - EPS = ( BASE**( 1-IT ) ) / 2 - ELSE - RND = ZERO - EPS = BASE**( 1-IT ) - END IF - PREC = EPS*BASE - EMIN = IMIN - EMAX = IMAX - SFMIN = RMIN - SMALL = ONE / RMAX - IF( SMALL.GE.SFMIN ) THEN -* -* Use SMALL plus a bit, to avoid the possibility of rounding -* causing overflow when computing 1/sfmin. -* - SFMIN = SMALL*( ONE+EPS ) - END IF - END IF -* - IF( LSAME( CMACH, 'E' ) ) THEN - RMACH = EPS - ELSE IF( LSAME( CMACH, 'S' ) ) THEN - RMACH = SFMIN - ELSE IF( LSAME( CMACH, 'B' ) ) THEN - RMACH = BASE - ELSE IF( LSAME( CMACH, 'P' ) ) THEN - RMACH = PREC - ELSE IF( LSAME( CMACH, 'N' ) ) THEN - RMACH = T - ELSE IF( LSAME( CMACH, 'R' ) ) THEN - RMACH = RND - ELSE IF( LSAME( CMACH, 'M' ) ) THEN - RMACH = EMIN - ELSE IF( LSAME( CMACH, 'U' ) ) THEN - RMACH = RMIN - ELSE IF( LSAME( CMACH, 'L' ) ) THEN - RMACH = EMAX - ELSE IF( LSAME( CMACH, 'O' ) ) THEN - RMACH = RMAX - END IF -* - QLAMCH = RMACH - RETURN -* -* End of QLAMCH -* - END -* -************************************************************************ -* - SUBROUTINE QLAMC1( BETA, T, RND, IEEE1 ) -* -* -- 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 .. - LOGICAL IEEE1, RND - INTEGER BETA, T -* .. -* -* Purpose -* ======= -* -* QLAMC1 determines the machine parameters given by BETA, T, RND, and -* IEEE1. -* -* Arguments -* ========= -* -* BETA (output) INTEGER -* The base of the machine. -* -* T (output) INTEGER -* The number of ( BETA ) digits in the mantissa. -* -* RND (output) LOGICAL -* Specifies whether proper rounding ( RND = .TRUE. ) or -* chopping ( RND = .FALSE. ) occurs in addition. This may not -* be a reliable guide to the way in which the machine performs -* its arithmetic. -* -* IEEE1 (output) LOGICAL -* Specifies whether rounding appears to be done in the IEEE -* 'round to nearest' style. -* -* Further Details -* =============== -* -* The routine is based on the routine ENVRON by Malcolm and -* incorporates suggestions by Gentleman and Marovich. See -* -* Malcolm M. A. (1972) Algorithms to reveal properties of -* floating-point arithmetic. Comms. of the ACM, 15, 949-951. -* -* Gentleman W. M. and Marovich S. B. (1974) More on algorithms -* that reveal properties of floating point arithmetic units. -* Comms. of the ACM, 17, 276-277. -* -* ===================================================================== -* -* .. Local Scalars .. - LOGICAL FIRST, LIEEE1, LRND - INTEGER LBETA, LT - REAL(16) A, B, C, F, ONE, QTR, SAVEC, T1, T2 -* .. -* .. External Functions .. - REAL(16) QLAMC3 - EXTERNAL QLAMC3 -* .. -* .. Save statement .. - SAVE FIRST, LIEEE1, LBETA, LRND, LT -* .. -* .. Data statements .. - DATA FIRST / .TRUE. / -* .. -* .. Executable Statements .. -* - IF( FIRST ) THEN - FIRST = .FALSE. - ONE = 1 -* -* LBETA, LIEEE1, LT and LRND are the local values of BETA, -* IEEE1, T and RND. -* -* Throughout this routine we use the function QLAMC3 to ensure -* that relevant values are stored and not held in registers, or -* are not affected by optimizers. -* -* Compute a = 2.0**m with the smallest positive integer m such -* that -* -* fl( a + 1.0 ) = a. -* - A = 1 - C = 1 -* -*+ WHILE( C.EQ.ONE )LOOP - 10 CONTINUE - IF( C.EQ.ONE ) THEN - A = 2*A - C = QLAMC3( A, ONE ) - C = QLAMC3( C, -A ) - GO TO 10 - END IF -*+ END WHILE -* -* Now compute b = 2.0**m with the smallest positive integer m -* such that -* -* fl( a + b ) .gt. a. -* - B = 1 - C = QLAMC3( A, B ) -* -*+ WHILE( C.EQ.A )LOOP - 20 CONTINUE - IF( C.EQ.A ) THEN - B = 2*B - C = QLAMC3( A, B ) - GO TO 20 - END IF -*+ END WHILE -* -* Now compute the base. a and c are neighbouring floating point -* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so -* their difference is beta. Adding 0.25 to c is to ensure that it -* is truncated to beta and not ( beta - 1 ). -* - QTR = ONE / 4 - SAVEC = C - C = QLAMC3( C, -A ) - LBETA = C + QTR -* -* Now determine whether rounding or chopping occurs, by adding a -* bit less than beta/2 and a bit more than beta/2 to a. -* - B = LBETA - F = QLAMC3( B / 2, -B / 100 ) - C = QLAMC3( F, A ) - IF( C.EQ.A ) THEN - LRND = .TRUE. - ELSE - LRND = .FALSE. - END IF - F = QLAMC3( B / 2, B / 100 ) - C = QLAMC3( F, A ) - IF( ( LRND ) .AND. ( C.EQ.A ) ) - $ LRND = .FALSE. -* -* Try and decide whether rounding is done in the IEEE 'round to -* nearest' style. B/2 is half a unit in the last place of the two -* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit -* zero, and SAVEC is odd. Thus adding B/2 to A should not change -* A, but adding B/2 to SAVEC should change SAVEC. -* - T1 = QLAMC3( B / 2, A ) - T2 = QLAMC3( B / 2, SAVEC ) - LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND -* -* Now find the mantissa, t. It should be the integer part of -* log to the base beta of a, however it is safer to determine t -* by powering. So we find t as the smallest positive integer for -* which -* -* fl( beta**t + 1.0 ) = 1.0. -* - LT = 0 - A = 1 - C = 1 -* -*+ WHILE( C.EQ.ONE )LOOP - 30 CONTINUE - IF( C.EQ.ONE ) THEN - LT = LT + 1 - A = A*LBETA - C = QLAMC3( A, ONE ) - C = QLAMC3( C, -A ) - GO TO 30 - END IF -*+ END WHILE -* - END IF -* - BETA = LBETA - T = LT - RND = LRND - IEEE1 = LIEEE1 - RETURN -* -* End of QLAMC1 -* - END -* -************************************************************************ -* - SUBROUTINE QLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) -* -* -- 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 .. - LOGICAL RND - INTEGER BETA, EMAX, EMIN, T - REAL(16) EPS, RMAX, RMIN -* .. -* -* Purpose -* ======= -* -* QLAMC2 determines the machine parameters specified in its argument -* list. -* -* Arguments -* ========= -* -* BETA (output) INTEGER -* The base of the machine. -* -* T (output) INTEGER -* The number of ( BETA ) digits in the mantissa. -* -* RND (output) LOGICAL -* Specifies whether proper rounding ( RND = .TRUE. ) or -* chopping ( RND = .FALSE. ) occurs in addition. This may not -* be a reliable guide to the way in which the machine performs -* its arithmetic. -* -* EPS (output) REAL(16) -* The smallest positive number such that -* -* fl( 1.0 - EPS ) .LT. 1.0, -* -* where fl denotes the computed value. -* -* EMIN (output) INTEGER -* The minimum exponent before (gradual) underflow occurs. -* -* RMIN (output) REAL(16) -* The smallest normalized number for the machine, given by -* BASE**( EMIN - 1 ), where BASE is the floating point value -* of BETA. -* -* EMAX (output) INTEGER -* The maximum exponent before overflow occurs. -* -* RMAX (output) REAL(16) -* The largest positive number for the machine, given by -* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point -* value of BETA. -* -* Further Details -* =============== -* -* The computation of EPS is based on a routine PARANOIA by -* W. Kahan of the University of California at Berkeley. -* -* ===================================================================== -* -* .. Local Scalars .. - LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND - INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, - $ NGNMIN, NGPMIN - REAL(16) A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, - $ SIXTH, SMALL, THIRD, TWO, ZERO -* .. -* .. External Functions .. - REAL(16) QLAMC3 - EXTERNAL QLAMC3 -* .. -* .. External Subroutines .. - EXTERNAL QLAMC1, QLAMC4, QLAMC5 -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, MAX, MIN -* .. -* .. Save statement .. - SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, - $ LRMIN, LT -* .. -* .. Data statements .. - DATA FIRST / .TRUE. / , IWARN / .FALSE. / -* .. -* .. Executable Statements .. -* - IF( FIRST ) THEN - FIRST = .FALSE. - ZERO = 0 - ONE = 1 - TWO = 2 -* -* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of -* BETA, T, RND, EPS, EMIN and RMIN. -* -* Throughout this routine we use the function QLAMC3 to ensure -* that relevant values are stored and not held in registers, or -* are not affected by optimizers. -* -* QLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. -* - CALL QLAMC1( LBETA, LT, LRND, LIEEE1 ) -* -* Start to find EPS. -* - B = LBETA - A = B**( -LT ) - LEPS = A -* -* Try some tricks to see whether or not this is the correct EPS. -* - B = TWO / 3 - HALF = ONE / 2 - SIXTH = QLAMC3( B, -HALF ) - THIRD = QLAMC3( SIXTH, SIXTH ) - B = QLAMC3( THIRD, -HALF ) - B = QLAMC3( B, SIXTH ) - B = ABS( B ) - IF( B.LT.LEPS ) - $ B = LEPS -* - LEPS = 1 -* -*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP - 10 CONTINUE - IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN - LEPS = B - C = QLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) - C = QLAMC3( HALF, -C ) - B = QLAMC3( HALF, C ) - C = QLAMC3( HALF, -B ) - B = QLAMC3( HALF, C ) - GO TO 10 - END IF -*+ END WHILE -* - IF( A.LT.LEPS ) - $ LEPS = A -* -* Computation of EPS complete. -* -* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). -* Keep dividing A by BETA until (gradual) underflow occurs. This -* is detected when we cannot recover the previous A. -* - RBASE = ONE / LBETA - SMALL = ONE - DO 20 I = 1, 3 - SMALL = QLAMC3( SMALL*RBASE, ZERO ) - 20 CONTINUE - A = QLAMC3( ONE, SMALL ) - CALL QLAMC4( NGPMIN, ONE, LBETA ) - CALL QLAMC4( NGNMIN, -ONE, LBETA ) - CALL QLAMC4( GPMIN, A, LBETA ) - CALL QLAMC4( GNMIN, -A, LBETA ) - IEEE = .FALSE. -* - IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN - IF( NGPMIN.EQ.GPMIN ) THEN - LEMIN = NGPMIN -* ( Non twos-complement machines, no gradual underflow; -* e.g., VAX ) - ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN - LEMIN = NGPMIN - 1 + LT - IEEE = .TRUE. -* ( Non twos-complement machines, with gradual underflow; -* e.g., IEEE standard followers ) - ELSE - LEMIN = MIN( NGPMIN, GPMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -* - ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN - IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN - LEMIN = MAX( NGPMIN, NGNMIN ) -* ( Twos-complement machines, no gradual underflow; -* e.g., CYBER 205 ) - ELSE - LEMIN = MIN( NGPMIN, NGNMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -* - ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. - $ ( GPMIN.EQ.GNMIN ) ) THEN - IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN - LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT -* ( Twos-complement machines with gradual underflow; -* no known machine ) - ELSE - LEMIN = MIN( NGPMIN, NGNMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -* - ELSE - LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -*** -* Comment out this if block if EMIN is ok - IF( IWARN ) THEN - FIRST = .TRUE. - WRITE( 6, FMT = 9999 )LEMIN - END IF -*** -* -* Assume IEEE arithmetic if we found denormalised numbers above, -* or if arithmetic seems to round in the IEEE style, determined -* in routine QLAMC1. A true IEEE machine should have both things -* true; however, faulty machines may have one or the other. -* - IEEE = IEEE .OR. LIEEE1 -* -* Compute RMIN by successive division by BETA. We could compute -* RMIN as BASE**( EMIN - 1 ), but some machines underflow during -* this computation. -* - LRMIN = 1 - DO 30 I = 1, 1 - LEMIN - LRMIN = QLAMC3( LRMIN*RBASE, ZERO ) - 30 CONTINUE -* -* Finally, call QLAMC5 to compute EMAX and RMAX. -* - CALL QLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) - END IF -* - BETA = LBETA - T = LT - RND = LRND - EPS = LEPS - EMIN = LEMIN - RMIN = LRMIN - EMAX = LEMAX - RMAX = LRMAX -* - RETURN -* - 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', - $ ' EMIN = ', I8, / - $ ' If, after inspection, the value EMIN looks', - $ ' acceptable please comment out ', - $ / ' the IF block as marked within the code of routine', - $ ' QLAMC2,', / ' otherwise supply EMIN explicitly.', / ) -* -* End of QLAMC2 -* - END -* -************************************************************************ -* - REAL(16) FUNCTION QLAMC3( A, B ) -* -* -- 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 .. - REAL(16) A, B -* .. -* -* Purpose -* ======= -* -* QLAMC3 is intended to force A and B to be stored prior to doing -* the addition of A and B , for use in situations where optimizers -* might hold one of these in a register. -* -* Arguments -* ========= -* -* A, B (input) REAL(16) -* The values A and B. -* -* ===================================================================== -* -* .. Executable Statements .. -* - QLAMC3 = A + B -* - RETURN -* -* End of QLAMC3 -* - END -* -************************************************************************ -* - SUBROUTINE QLAMC4( EMIN, START, BASE ) -* -* -- 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 BASE, EMIN - REAL(16) START -* .. -* -* Purpose -* ======= -* -* QLAMC4 is a service routine for QLAMC2. -* -* Arguments -* ========= -* -* EMIN (output) EMIN -* The minimum exponent before (gradual) underflow, computed by -* setting A = START and dividing by BASE until the previous A -* can not be recovered. -* -* START (input) REAL(16) -* The starting point for determining EMIN. -* -* BASE (input) INTEGER -* The base of the machine. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I - REAL(16) A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO -* .. -* .. External Functions .. - REAL(16) QLAMC3 - EXTERNAL QLAMC3 -* .. -* .. Executable Statements .. -* - A = START - ONE = 1 - RBASE = ONE / BASE - ZERO = 0 - EMIN = 1 - B1 = QLAMC3( A*RBASE, ZERO ) - C1 = A - C2 = A - D1 = A - D2 = A -*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. -* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP - 10 CONTINUE - IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. - $ ( D2.EQ.A ) ) THEN - EMIN = EMIN - 1 - A = B1 - B1 = QLAMC3( A / BASE, ZERO ) - C1 = QLAMC3( B1*BASE, ZERO ) - D1 = ZERO - DO 20 I = 1, BASE - D1 = D1 + B1 - 20 CONTINUE - B2 = QLAMC3( A*RBASE, ZERO ) - C2 = QLAMC3( B2 / RBASE, ZERO ) - D2 = ZERO - DO 30 I = 1, BASE - D2 = D2 + B2 - 30 CONTINUE - GO TO 10 - END IF -*+ END WHILE -* - RETURN -* -* End of QLAMC4 -* - END -* -************************************************************************ -* - SUBROUTINE QLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) -* -* -- 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 .. - LOGICAL IEEE - INTEGER BETA, EMAX, EMIN, P - REAL(16) RMAX -* .. -* -* Purpose -* ======= -* -* QLAMC5 attempts to compute RMAX, the largest machine floating-point -* number, without overflow. It assumes that EMAX + abs(EMIN) sum -* approximately to a power of 2. It will fail on machines where this -* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, -* EMAX = 28718). It will also fail if the value supplied for EMIN is -* too large (i.e. too close to zero), probably with overflow. -* -* Arguments -* ========= -* -* BETA (input) INTEGER -* The base of floating-point arithmetic. -* -* P (input) INTEGER -* The number of base BETA digits in the mantissa of a -* floating-point value. -* -* EMIN (input) INTEGER -* The minimum exponent before (gradual) underflow. -* -* IEEE (input) LOGICAL -* A logical flag specifying whether or not the arithmetic -* system is thought to comply with the IEEE standard. -* -* EMAX (output) INTEGER -* The largest exponent before overflow -* -* RMAX (output) REAL(16) -* The largest machine floating-point number. -* -* ===================================================================== -* -* .. Parameters .. - REAL(16) ZERO, ONE - PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) -* .. -* .. Local Scalars .. - INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP - REAL(16) OLDY, RECBAS, Y, Z -* .. -* .. External Functions .. - REAL(16) QLAMC3 - EXTERNAL QLAMC3 -* .. -* .. Intrinsic Functions .. - INTRINSIC MOD -* .. -* .. Executable Statements .. -* -* First compute LEXP and UEXP, two powers of 2 that bound -* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum -* approximately to the bound that is closest to abs(EMIN). -* (EMAX is the exponent of the required number RMAX). -* - LEXP = 1 - EXBITS = 1 - 10 CONTINUE - TRY = LEXP*2 - IF( TRY.LE.( -EMIN ) ) THEN - LEXP = TRY - EXBITS = EXBITS + 1 - GO TO 10 - END IF - IF( LEXP.EQ.-EMIN ) THEN - UEXP = LEXP - ELSE - UEXP = TRY - EXBITS = EXBITS + 1 - END IF -* -* Now -LEXP is less than or equal to EMIN, and -UEXP is greater -* than or equal to EMIN. EXBITS is the number of bits needed to -* store the exponent. -* - IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN - EXPSUM = 2*LEXP - ELSE - EXPSUM = 2*UEXP - END IF -* -* EXPSUM is the exponent range, approximately equal to -* EMAX - EMIN + 1 . -* - EMAX = EXPSUM + EMIN - 1 - NBITS = 1 + EXBITS + P -* -* NBITS is the total number of bits needed to store a -* floating-point number. -* - IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN -* -* Either there are an odd number of bits used to store a -* floating-point number, which is unlikely, or some bits are -* not used in the representation of numbers, which is possible, -* (e.g. Cray machines) or the mantissa has an implicit bit, -* (e.g. IEEE machines, Dec Vax machines), which is perhaps the -* most likely. We have to assume the last alternative. -* If this is true, then we need to reduce EMAX by one because -* there must be some way of representing zero in an implicit-bit -* system. On machines like Cray, we are reducing EMAX by one -* unnecessarily. -* - EMAX = EMAX - 1 - END IF -* - IF( IEEE ) THEN -* -* Assume we are on an IEEE machine which reserves one exponent -* for infinity and NaN. -* - EMAX = EMAX - 1 - END IF -* -* Now create RMAX, the largest machine number, which should -* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . -* -* First compute 1.0 - BETA**(-P), being careful that the -* result is less than 1.0 . -* - RECBAS = ONE / BETA - Z = BETA - ONE - Y = ZERO - DO 20 I = 1, P - Z = Z*RECBAS - IF( Y.LT.ONE ) - $ OLDY = Y - Y = QLAMC3( Y, Z ) - 20 CONTINUE - IF( Y.GE.ONE ) - $ Y = OLDY -* -* Now multiply by BETA**EMAX to get RMAX. -* - DO 30 I = 1, EMAX - Y = QLAMC3( Y*BETA, ZERO ) - 30 CONTINUE -* - RMAX = Y - RETURN -* -* End of QLAMC5 -* - END diff --git a/mtx/lapack_quad/qlaswp.f b/mtx/lapack_quad/qlaswp.f deleted file mode 100644 index 2567521ba..000000000 --- a/mtx/lapack_quad/qlaswp.f +++ /dev/null @@ -1,119 +0,0 @@ - SUBROUTINE QLASWP( N, A, LDA, K1, K2, IPIV, INCX ) -* -* -- LAPACK auxiliary routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. Scalar Arguments .. - INTEGER INCX, K1, K2, LDA, N -* .. -* .. Array Arguments .. - INTEGER IPIV( * ) - REAL(16) A( LDA, * ) -* .. -* -* Purpose -* ======= -* -* QLASWP performs a series of row interchanges on the matrix A. -* One row interchange is initiated for each of rows K1 through K2 of A. -* -* Arguments -* ========= -* -* N (input) INTEGER -* The number of columns of the matrix A. -* -* A (input/output) REAL(16) array, dimension (LDA,N) -* On entry, the matrix of column dimension N to which the row -* interchanges will be applied. -* On exit, the permuted matrix. -* -* LDA (input) INTEGER -* The leading dimension of the array A. -* -* K1 (input) INTEGER -* The first element of IPIV for which a row interchange will -* be done. -* -* K2 (input) INTEGER -* The last element of IPIV for which a row interchange will -* be done. -* -* IPIV (input) INTEGER array, dimension (K2*abs(INCX)) -* The vector of pivot indices. Only the elements in positions -* K1 through K2 of IPIV are accessed. -* IPIV(K) = L implies rows K and L are to be interchanged. -* -* INCX (input) INTEGER -* The increment between successive values of IPIV. If IPIV -* is negative, the pivots are applied in reverse order. -* -* Further Details -* =============== -* -* Modified by -* R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I, I1, I2, INC, IP, IX, IX0, J, K, N32 - REAL(16) TEMP -* .. -* .. Executable Statements .. -* -* Interchange row I with row IPIV(I) for each of rows K1 through K2. -* - IF( INCX.GT.0 ) THEN - IX0 = K1 - I1 = K1 - I2 = K2 - INC = 1 - ELSE IF( INCX.LT.0 ) THEN - IX0 = 1 + ( 1-K2 )*INCX - I1 = K2 - I2 = K1 - INC = -1 - ELSE - RETURN - END IF -* - N32 = ( N / 32 )*32 - IF( N32.NE.0 ) THEN - DO 30 J = 1, N32, 32 - IX = IX0 - DO 20 I = I1, I2, INC - IP = IPIV( IX ) - IF( IP.NE.I ) THEN - DO 10 K = J, J + 31 - TEMP = A( I, K ) - A( I, K ) = A( IP, K ) - A( IP, K ) = TEMP - 10 CONTINUE - END IF - IX = IX + INCX - 20 CONTINUE - 30 CONTINUE - END IF - IF( N32.NE.N ) THEN - N32 = N32 + 1 - IX = IX0 - DO 50 I = I1, I2, INC - IP = IPIV( IX ) - IF( IP.NE.I ) THEN - DO 40 K = N32, N - TEMP = A( I, K ) - A( I, K ) = A( IP, K ) - A( IP, K ) = TEMP - 40 CONTINUE - END IF - IX = IX + INCX - 50 CONTINUE - END IF -* - RETURN -* -* End of QLASWP -* - END diff --git a/mtx/lapack_quad/qscal.f b/mtx/lapack_quad/qscal.f deleted file mode 100644 index 4b6f1e1f4..000000000 --- a/mtx/lapack_quad/qscal.f +++ /dev/null @@ -1,57 +0,0 @@ - SUBROUTINE QSCAL(N,DA,DX,INCX) -* .. Scalar Arguments .. - REAL(16) DA - INTEGER INCX,N -* .. -* .. Array Arguments .. - REAL(16) DX(*) -* .. -* -* Purpose -* ======= -** -* scales a vector by a constant. -* uses unrolled loops for increment equal to one. -* jack dongarra, linpack, 3/11/78. -* modified 3/93 to return if incx .le. 0. -* modified 12/3/93, array(1) declarations changed to array(*) -* -* -* .. Local Scalars .. - INTEGER I,M,MP1,NINCX -* .. -* .. Intrinsic Functions .. - INTRINSIC MOD -* .. - IF (N.LE.0 .OR. INCX.LE.0) RETURN - IF (INCX.EQ.1) GO TO 20 -* -* code for increment not equal to 1 -* - NINCX = N*INCX - DO 10 I = 1,NINCX,INCX - DX(I) = DA*DX(I) - 10 CONTINUE - RETURN -* -* code for increment equal to 1 -* -* -* clean-up loop -* - 20 M = MOD(N,5) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DX(I) = DA*DX(I) - 30 CONTINUE - IF (N.LT.5) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,5 - DX(I) = DA*DX(I) - DX(I+1) = DA*DX(I+1) - DX(I+2) = DA*DX(I+2) - DX(I+3) = DA*DX(I+3) - DX(I+4) = DA*DX(I+4) - 50 CONTINUE - RETURN - END diff --git a/mtx/lapack_quad/qswap.f b/mtx/lapack_quad/qswap.f deleted file mode 100644 index 556bb6cc8..000000000 --- a/mtx/lapack_quad/qswap.f +++ /dev/null @@ -1,70 +0,0 @@ - SUBROUTINE QSWAP(N,DX,INCX,DY,INCY) -* .. Scalar Arguments .. - INTEGER INCX,INCY,N -* .. -* .. Array Arguments .. - REAL(16) DX(*),DY(*) -* .. -* -* Purpose -* ======= -* -* interchanges two vectors. -* uses unrolled loops for increments equal one. -* jack dongarra, linpack, 3/11/78. -* modified 12/3/93, array(1) declarations changed to array(*) -* -* -* .. Local Scalars .. - REAL(16) DTEMP - INTEGER I,IX,IY,M,MP1 -* .. -* .. Intrinsic Functions .. - INTRINSIC MOD -* .. - 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 - DTEMP = DX(IX) - DX(IX) = DY(IY) - DY(IY) = DTEMP - IX = IX + INCX - IY = IY + INCY - 10 CONTINUE - RETURN -* -* code for both increments equal to 1 -* -* -* clean-up loop -* - 20 M = MOD(N,3) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DTEMP = DX(I) - DX(I) = DY(I) - DY(I) = DTEMP - 30 CONTINUE - IF (N.LT.3) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,3 - DTEMP = DX(I) - DX(I) = DY(I) - DY(I) = DTEMP - DTEMP = DX(I+1) - DX(I+1) = DY(I+1) - DY(I+1) = DTEMP - DTEMP = DX(I+2) - DX(I+2) = DY(I+2) - DY(I+2) = DTEMP - 50 CONTINUE - RETURN - END diff --git a/mtx/lapack_quad/qtrsm.f b/mtx/lapack_quad/qtrsm.f deleted file mode 100644 index 44c0fb5cb..000000000 --- a/mtx/lapack_quad/qtrsm.f +++ /dev/null @@ -1,373 +0,0 @@ - SUBROUTINE QTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) -* .. Scalar Arguments .. - REAL(16) ALPHA - INTEGER LDA,LDB,M,N - CHARACTER DIAG,SIDE,TRANSA,UPLO -* .. -* .. Array Arguments .. - REAL(16) A(LDA,*),B(LDB,*) -* .. -* -* Purpose -* ======= -* -* QTRSM solves one of the matrix equations -* -* op( A )*X = alpha*B, or X*op( A ) = alpha*B, -* -* where alpha is a scalar, X and B are m by n matrices, A is a unit, or -* non-unit, upper or lower triangular matrix and op( A ) is one of -* -* op( A ) = A or op( A ) = A'. -* -* The matrix X is overwritten on B. -* -* Arguments -* ========== -* -* SIDE - CHARACTER*1. -* On entry, SIDE specifies whether op( A ) appears on the left -* or right of X as follows: -* -* SIDE = 'L' or 'l' op( A )*X = alpha*B. -* -* SIDE = 'R' or 'r' X*op( A ) = alpha*B. -* -* Unchanged on exit. -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix A 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. -* -* TRANSA - CHARACTER*1. -* On entry, TRANSA specifies the form of op( A ) to be used in -* the matrix multiplication as follows: -* -* TRANSA = 'N' or 'n' op( A ) = A. -* -* TRANSA = 'T' or 't' op( A ) = A'. -* -* TRANSA = 'C' or 'c' op( A ) = A'. -* -* 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. -* -* M - INTEGER. -* On entry, M specifies the number of rows of B. M must be at -* least zero. -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the number of columns of B. N must be -* at least zero. -* Unchanged on exit. -* -* ALPHA - REAL(16). -* On entry, ALPHA specifies the scalar alpha. When alpha is -* zero then A is not referenced and B need not be set before -* entry. -* Unchanged on exit. -* -* A - REAL(16) array of DIMENSION ( LDA, k ), where k is m -* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. -* Before entry with UPLO = 'U' or 'u', the leading k by k -* upper triangular part of the array A must contain the upper -* triangular matrix and the strictly lower triangular part of -* A is not referenced. -* Before entry with UPLO = 'L' or 'l', the leading k by k -* lower triangular part of the array A must contain the lower -* triangular matrix and the strictly upper triangular part of -* A is not referenced. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced either, 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. When SIDE = 'L' or 'l' then -* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' -* then LDA must be at least max( 1, n ). -* Unchanged on exit. -* -* B - REAL(16) array of DIMENSION ( LDB, n ). -* Before entry, the leading m by n part of the array B must -* contain the right-hand side matrix B, and on exit is -* overwritten by the solution matrix X. -* -* LDB - INTEGER. -* On entry, LDB specifies the first dimension of B as declared -* in the calling (sub) program. LDB must be at least -* max( 1, m ). -* Unchanged on exit. -* -* -* Level 3 Blas routine. -* -* -* -- Written on 8-February-1989. -* Jack Dongarra, Argonne National Laboratory. -* Iain Duff, AERE Harwell. -* Jeremy Du Croz, Numerical Algorithms Group Ltd. -* Sven Hammarling, Numerical Algorithms Group Ltd. -* -* -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* .. Local Scalars .. - REAL(16) TEMP - INTEGER I,INFO,J,K,NROWA - LOGICAL LSIDE,NOUNIT,UPPER -* .. -* .. Parameters .. - REAL(16) ONE,ZERO - PARAMETER (ONE=1.0_16,ZERO=0.0_16) -* .. -* -* Test the input parameters. -* - LSIDE = LSAME(SIDE,'L') - IF (LSIDE) THEN - NROWA = M - ELSE - NROWA = N - END IF - NOUNIT = LSAME(DIAG,'N') - UPPER = LSAME(UPLO,'U') -* - INFO = 0 - IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN - INFO = 1 - ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN - INFO = 2 - ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND. - + (.NOT.LSAME(TRANSA,'T')) .AND. - + (.NOT.LSAME(TRANSA,'C'))) THEN - INFO = 3 - ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN - INFO = 4 - ELSE IF (M.LT.0) THEN - INFO = 5 - ELSE IF (N.LT.0) THEN - INFO = 6 - ELSE IF (LDA.LT.MAX(1,NROWA)) THEN - INFO = 9 - ELSE IF (LDB.LT.MAX(1,M)) THEN - INFO = 11 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('QTRSM ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* -* And when alpha.eq.zero. -* - IF (ALPHA.EQ.ZERO) THEN - DO 20 J = 1,N - DO 10 I = 1,M - B(I,J) = ZERO - 10 CONTINUE - 20 CONTINUE - RETURN - END IF -* -* Start the operations. -* - IF (LSIDE) THEN - IF (LSAME(TRANSA,'N')) THEN -* -* Form B := alpha*inv( A )*B. -* - IF (UPPER) THEN - DO 60 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 30 I = 1,M - B(I,J) = ALPHA*B(I,J) - 30 CONTINUE - END IF - DO 50 K = M,1,-1 - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 40 I = 1,K - 1 - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 40 CONTINUE - END IF - 50 CONTINUE - 60 CONTINUE - ELSE - DO 100 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 70 I = 1,M - B(I,J) = ALPHA*B(I,J) - 70 CONTINUE - END IF - DO 90 K = 1,M - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 80 I = K + 1,M - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 80 CONTINUE - END IF - 90 CONTINUE - 100 CONTINUE - END IF - ELSE -* -* Form B := alpha*inv( A' )*B. -* - IF (UPPER) THEN - DO 130 J = 1,N - DO 120 I = 1,M - TEMP = ALPHA*B(I,J) - DO 110 K = 1,I - 1 - TEMP = TEMP - A(K,I)*B(K,J) - 110 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(I,I) - B(I,J) = TEMP - 120 CONTINUE - 130 CONTINUE - ELSE - DO 160 J = 1,N - DO 150 I = M,1,-1 - TEMP = ALPHA*B(I,J) - DO 140 K = I + 1,M - TEMP = TEMP - A(K,I)*B(K,J) - 140 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(I,I) - B(I,J) = TEMP - 150 CONTINUE - 160 CONTINUE - END IF - END IF - ELSE - IF (LSAME(TRANSA,'N')) THEN -* -* Form B := alpha*B*inv( A ). -* - IF (UPPER) THEN - DO 210 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 170 I = 1,M - B(I,J) = ALPHA*B(I,J) - 170 CONTINUE - END IF - DO 190 K = 1,J - 1 - IF (A(K,J).NE.ZERO) THEN - DO 180 I = 1,M - B(I,J) = B(I,J) - A(K,J)*B(I,K) - 180 CONTINUE - END IF - 190 CONTINUE - IF (NOUNIT) THEN - TEMP = ONE/A(J,J) - DO 200 I = 1,M - B(I,J) = TEMP*B(I,J) - 200 CONTINUE - END IF - 210 CONTINUE - ELSE - DO 260 J = N,1,-1 - IF (ALPHA.NE.ONE) THEN - DO 220 I = 1,M - B(I,J) = ALPHA*B(I,J) - 220 CONTINUE - END IF - DO 240 K = J + 1,N - IF (A(K,J).NE.ZERO) THEN - DO 230 I = 1,M - B(I,J) = B(I,J) - A(K,J)*B(I,K) - 230 CONTINUE - END IF - 240 CONTINUE - IF (NOUNIT) THEN - TEMP = ONE/A(J,J) - DO 250 I = 1,M - B(I,J) = TEMP*B(I,J) - 250 CONTINUE - END IF - 260 CONTINUE - END IF - ELSE -* -* Form B := alpha*B*inv( A' ). -* - IF (UPPER) THEN - DO 310 K = N,1,-1 - IF (NOUNIT) THEN - TEMP = ONE/A(K,K) - DO 270 I = 1,M - B(I,K) = TEMP*B(I,K) - 270 CONTINUE - END IF - DO 290 J = 1,K - 1 - IF (A(J,K).NE.ZERO) THEN - TEMP = A(J,K) - DO 280 I = 1,M - B(I,J) = B(I,J) - TEMP*B(I,K) - 280 CONTINUE - END IF - 290 CONTINUE - IF (ALPHA.NE.ONE) THEN - DO 300 I = 1,M - B(I,K) = ALPHA*B(I,K) - 300 CONTINUE - END IF - 310 CONTINUE - ELSE - DO 360 K = 1,N - IF (NOUNIT) THEN - TEMP = ONE/A(K,K) - DO 320 I = 1,M - B(I,K) = TEMP*B(I,K) - 320 CONTINUE - END IF - DO 340 J = K + 1,N - IF (A(J,K).NE.ZERO) THEN - TEMP = A(J,K) - DO 330 I = 1,M - B(I,J) = B(I,J) - TEMP*B(I,K) - 330 CONTINUE - END IF - 340 CONTINUE - IF (ALPHA.NE.ONE) THEN - DO 350 I = 1,M - B(I,K) = ALPHA*B(I,K) - 350 CONTINUE - END IF - 360 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of QTRSM . -* - END diff --git a/mtx/lapack_quad/qtrsv.f b/mtx/lapack_quad/qtrsv.f deleted file mode 100644 index 2d3b13c60..000000000 --- a/mtx/lapack_quad/qtrsv.f +++ /dev/null @@ -1,338 +0,0 @@ -*> \brief \b DTRSV -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -* Definition: -* =========== -* -* SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) -* -* .. Scalar Arguments .. -* INTEGER INCX,LDA,N -* CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. -* REAL(16) A(LDA,*),X(*) -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> DTRSV solves one of the systems of equations -*> -*> A*x = b, or A**T*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 matrix. -*> -*> No test for singularity or near-singularity is included in this -*> routine. Such tests must be performed before calling this routine. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] UPLO -*> \verbatim -*> UPLO is 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. -*> \endverbatim -*> -*> \param[in] TRANS -*> \verbatim -*> TRANS is 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**T*x = b. -*> -*> TRANS = 'C' or 'c' A**T*x = b. -*> \endverbatim -*> -*> \param[in] DIAG -*> \verbatim -*> DIAG is 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. -*> \endverbatim -*> -*> \param[in] N -*> \verbatim -*> N is INTEGER -*> On entry, N specifies the order of the matrix A. -*> N must be at least zero. -*> \endverbatim -*> -*> \param[in] A -*> \verbatim -*> A is REAL(16) array of DIMENSION ( LDA, n ). -*> Before entry with UPLO = 'U' or 'u', the leading n by n -*> upper triangular part of the array A must contain the upper -*> triangular matrix and the strictly lower triangular part of -*> A is not referenced. -*> Before entry with UPLO = 'L' or 'l', the leading n by n -*> lower triangular part of the array A must contain the lower -*> triangular matrix and the strictly upper triangular part of -*> A is not referenced. -*> Note that when DIAG = 'U' or 'u', the diagonal elements of -*> A are not referenced either, but are assumed to be unity. -*> \endverbatim -*> -*> \param[in] LDA -*> \verbatim -*> LDA is INTEGER -*> On entry, LDA specifies the first dimension of A as declared -*> in the calling (sub) program. LDA must be at least -*> max( 1, n ). -*> \endverbatim -*> -*> \param[in,out] X -*> \verbatim -*> X is REAL(16) 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. -*> \endverbatim -*> -*> \param[in] INCX -*> \verbatim -*> INCX is INTEGER -*> On entry, INCX specifies the increment for the elements of -*> X. INCX must not be zero. -*> -*> 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. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \date November 2011 -* -*> \ingroup double_blas_level1 -* -* ===================================================================== - SUBROUTINE QTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) -* -* -- Reference BLAS level1 routine (version 3.4.0) -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2011 -* -* .. Scalar Arguments .. - INTEGER INCX,LDA,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - REAL(16) A(LDA,*),X(*) -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL(16) ZERO - PARAMETER (ZERO=0.0D+0) -* .. -* .. Local Scalars .. - REAL(16) TEMP - INTEGER I,INFO,IX,J,JX,KX - LOGICAL NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC MAX -* .. -* -* 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 (LDA.LT.MAX(1,N)) THEN - INFO = 6 - ELSE IF (INCX.EQ.0) THEN - INFO = 8 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('DTRSV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - 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 sequentially with one pass through A. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x := inv( A )*x. -* - IF (LSAME(UPLO,'U')) THEN - IF (INCX.EQ.1) THEN - DO 20 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/A(J,J) - TEMP = X(J) - DO 10 I = J - 1,1,-1 - X(I) = X(I) - TEMP*A(I,J) - 10 CONTINUE - END IF - 20 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 40 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/A(J,J) - TEMP = X(JX) - IX = JX - DO 30 I = J - 1,1,-1 - IX = IX - INCX - X(IX) = X(IX) - TEMP*A(I,J) - 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 - IF (NOUNIT) X(J) = X(J)/A(J,J) - TEMP = X(J) - DO 50 I = J + 1,N - X(I) = X(I) - TEMP*A(I,J) - 50 CONTINUE - END IF - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/A(J,J) - TEMP = X(JX) - IX = JX - DO 70 I = J + 1,N - IX = IX + INCX - X(IX) = X(IX) - TEMP*A(I,J) - 70 CONTINUE - END IF - JX = JX + INCX - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := inv( A**T )*x. -* - IF (LSAME(UPLO,'U')) THEN - IF (INCX.EQ.1) THEN - DO 100 J = 1,N - TEMP = X(J) - DO 90 I = 1,J - 1 - TEMP = TEMP - A(I,J)*X(I) - 90 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(J,J) - X(J) = TEMP - 100 CONTINUE - ELSE - JX = KX - DO 120 J = 1,N - TEMP = X(JX) - IX = KX - DO 110 I = 1,J - 1 - TEMP = TEMP - A(I,J)*X(IX) - IX = IX + INCX - 110 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(J,J) - X(JX) = TEMP - JX = JX + INCX - 120 CONTINUE - END IF - ELSE - IF (INCX.EQ.1) THEN - DO 140 J = N,1,-1 - TEMP = X(J) - DO 130 I = N,J + 1,-1 - TEMP = TEMP - A(I,J)*X(I) - 130 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(J,J) - X(J) = TEMP - 140 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 160 J = N,1,-1 - TEMP = X(JX) - IX = KX - DO 150 I = N,J + 1,-1 - TEMP = TEMP - A(I,J)*X(IX) - IX = IX - INCX - 150 CONTINUE - IF (NOUNIT) TEMP = TEMP/A(J,J) - X(JX) = TEMP - JX = JX - INCX - 160 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of DTRSV . -* - END diff --git a/mtx/make/makefile_base b/mtx/make/makefile_base index a288df3cd..f3c336a71 100644 --- a/mtx/make/makefile_base +++ b/mtx/make/makefile_base @@ -183,24 +183,6 @@ ifeq ($(WHICH_LAPACK),USE_SRCS) dlarfx.f endif -LAPACK_QUAD_SRCS = \ - qgttrf.f \ - qgttrs.f \ - qgtts2.f \ - qgetrf.f \ - qgetrs.f \ - qgemv.f \ - qgemm.f \ - qgetf2.f \ - qswap.f \ - iqamax.f \ - qtrsm.f \ - qtrsv.f \ - qger.f \ - qscal.f \ - qlamch.f \ - qlaswp.f - ifeq ($(WHICH_BLAS),USE_SRCS) BLAS_SRCS = \ dgemm.f \ @@ -236,23 +218,13 @@ endif MTX_SRCS = \ mtx_def.f \ - my_lapack95_dble.f90 \ - my_lapack95_quad.f90 \ + my_lapack95.f90 \ mtx_support.f90 \ DGBSVX_wrapper.f90 \ pre_conditioners.f90 \ bcyclic.f90 \ mtx_lib.f90 \ -ifneq ($(MAKECMDGOALS),clean) - $(shell $(CPP) -DDBLE $(MOD_PRIVATE_DIR)/my_lapack95.F90 > my_lapack95_dble.f90) - $(shell $(CPP) $(MOD_PRIVATE_DIR)/my_lapack95.F90 > my_lapack95_quad.f90) - # Alter the timestamps so make wont force a recompile just because we copied the files, - # instead we only recompile if the original .F90 file was changed - $(shell touch -r $(MOD_PRIVATE_DIR)/my_lapack95.F90 my_lapack95_dble.f90) - $(shell touch -r $(MOD_PRIVATE_DIR)/my_lapack95.F90 my_lapack95_quad.f90) -endif - ################################################################# # # TARGETS @@ -266,7 +238,7 @@ ifeq ($(WHICH_LAPACK),USE_SRCS) BLAS_LIB = libmesablas.$(LIB_SUFFIX) endif -MTX_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(MTX_SRCS) $(LAPACK_QUAD_SRCS))) +MTX_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(MTX_SRCS))) LAPACK_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(LAPACK_SRCS))) BLAS_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(BLAS_SRCS))) @@ -336,8 +308,6 @@ COMPILE_CMD = $(COMPILE) $(LAPACK_OBJS) : COMPILE_CMD = $(COMPILE_XTRA) -std=gnu -w -$(LAPACK_QUAD_OBJS) : COMPILE_CMD = $(COMPILE_XTRA) -w - $(filter-out dlamch.o,$(BLAS_OBJS)) : COMPILE_CMD = $(COMPILE_XTRA) -std=gnu -w # must turn off optimization for dlamch or can get infinite loop!!! @@ -368,7 +338,7 @@ endif # # DEPENDENCIES -SRC_PATH = $(MOD_PUBLIC_DIR):$(MOD_PRIVATE_DIR):../blas_src:../lapack_src:../lapack_quad +SRC_PATH = $(MOD_PUBLIC_DIR):$(MOD_PRIVATE_DIR):../blas_src:../lapack_src vpath %.f $(SRC_PATH) vpath %.f90 $(SRC_PATH) diff --git a/mtx/private/bcyclic.f90 b/mtx/private/bcyclic.f90 index 3d93f4bfa..d0b78dcca 100644 --- a/mtx/private/bcyclic.f90 +++ b/mtx/private/bcyclic.f90 @@ -23,7 +23,7 @@ module bcyclic use const_def, only: dp - use my_lapack95_dble + use my_lapack95 use utils_lib, only: set_nan, mesa_error implicit none diff --git a/mtx/private/mtx_support.f90 b/mtx/private/mtx_support.f90 index 41fdf71a0..3edf6f51c 100644 --- a/mtx/private/mtx_support.f90 +++ b/mtx/private/mtx_support.f90 @@ -666,7 +666,7 @@ end subroutine do_col_sparse_0_based_to_dense subroutine do_block_dble_mv(nvar, nz, lblk, dblk, ublk, b, prod) - use my_lapack95_dble, only: my_gemv_p1 + use my_lapack95, only: my_gemv_p1 integer, intent(in) :: nvar, nz real(dp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz) real(dp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz) diff --git a/mtx/private/my_lapack95.F90 b/mtx/private/my_lapack95.f90 similarity index 99% rename from mtx/private/my_lapack95.F90 rename to mtx/private/my_lapack95.f90 index f9ece650d..7401faa05 100644 --- a/mtx/private/my_lapack95.F90 +++ b/mtx/private/my_lapack95.f90 @@ -23,19 +23,11 @@ ! ! *********************************************************************** -#ifdef DBLE - module my_lapack95_dble + module my_lapack95 use const_def, only: dp use utils_lib, only: mesa_error implicit none integer, parameter :: fltp = dp -#else - module my_lapack95_quad - use const_def, only: qp - use utils_lib, only: mesa_error - implicit none - integer, parameter :: fltp = qp -#endif @@ -1432,9 +1424,4 @@ subroutine my_laswp_dbg( n, a, lda, k1, k2, ipiv, incx ) end if end subroutine my_laswp_dbg - -#ifdef DBLE - end module my_lapack95_dble -#else - end module my_lapack95_quad -#endif + end module my_lapack95 diff --git a/mtx/public/mtx_lapack95.dek b/mtx/public/mtx_lapack95.dek index f49b4ebd5..208c6d0f6 100644 --- a/mtx/public/mtx_lapack95.dek +++ b/mtx/public/mtx_lapack95.dek @@ -1,5 +1,5 @@ subroutine my_getf2_95(m, a, lda, ipiv, sfmin, info) - use my_lapack95_dble, only: my_getf2 + use my_lapack95, only: my_getf2 integer :: info, lda, m integer :: ipiv(:) real(dp) :: a(:,:) @@ -8,11 +8,9 @@ end subroutine my_getf2_95 subroutine my_getrs_95( n, nrhs, a, lda, ipiv, b, ldb, info ) - use my_lapack95_dble, only: my_getrs + use my_lapack95, only: my_getrs integer :: info, lda, ldb, n, nrhs integer, pointer :: ipiv(:) real(dp), pointer :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * ) call my_getrs( n, nrhs, a, lda, ipiv, b, ldb, info ) end subroutine my_getrs_95 - - \ No newline at end of file diff --git a/mtx/test/make/makefile_base b/mtx/test/make/makefile_base index 98e466243..356e90779 100644 --- a/mtx/test/make/makefile_base +++ b/mtx/test/make/makefile_base @@ -10,7 +10,7 @@ include $(MESA_DIR)/utils/makefile_header SRCS = \ test_mtx_support.f90 test_block_tri_dble.f90 test_block_tri_quad.f90 \ - test_square.f90 test_square_quad.f90 test_mtx.f90 + test_square.f90 test_mtx.f90 ################################################################# # diff --git a/mtx/test/src/test_mtx.f90 b/mtx/test/src/test_mtx.f90 index c13052af5..66173b251 100644 --- a/mtx/test/src/test_mtx.f90 +++ b/mtx/test/src/test_mtx.f90 @@ -27,7 +27,6 @@ program test_mtx use test_mtx_support use test_square - use test_square_quad use test_block_tri_dble, only: do_test_block_tri_dble use test_block_tri_quad, only: do_test_block_tri_quad @@ -49,7 +48,6 @@ program test_mtx call math_init() call do_test_square - call do_test_square_quad call do_test_block_tri_dble call do_test_block_tri_quad call test_format_conversion diff --git a/mtx/test/src/test_mtx_support.f90 b/mtx/test/src/test_mtx_support.f90 index cec81c1a7..525279c58 100644 --- a/mtx/test/src/test_mtx_support.f90 +++ b/mtx/test/src/test_mtx_support.f90 @@ -106,44 +106,4 @@ subroutine test_format_conversion end subroutine test_format_conversion - subroutine test_quad_tridiag - integer, parameter :: n = 5 - real(16), dimension(n) :: DL, D, DU, DU2, B - integer, dimension(n) :: ip - integer :: ierr, i - - write (*, *) 'test_quad_tridiag' - - DL = -2 ! subdiagonal - D = 3 ! diagonal - DU = -1 ! superdiagonal - - do i = 1, n - b(i) = i - 1 - end do - - ierr = 0 - ! factor - call qgttrf(n, DL, D, DU, DU2, ip, ierr) - if (ierr /= 0) then - write (*, *) 'failed in factoring' - call mesa_error(__FILE__, __LINE__) - end if - - ierr = 0 - ! solve - call qgttrs('N', n, 1, DL, D, DU, DU2, ip, B, n, ierr) - if (ierr /= 0) then - write (*, *) 'failed in solving' - call mesa_error(__FILE__, __LINE__) - end if - - do i = 1, n - write (*, *) i, b(i) - end do - write (*, *) - write (*, *) - - end subroutine test_quad_tridiag - end module test_mtx_support diff --git a/mtx/test/src/test_square_quad.f90 b/mtx/test/src/test_square_quad.f90 deleted file mode 100644 index 505a27b69..000000000 --- a/mtx/test/src/test_square_quad.f90 +++ /dev/null @@ -1,175 +0,0 @@ -! *********************************************************************** -! -! Copyright (C) 2011 Bill Paxton -! -! This file is part of MESA. -! -! MESA is free software; you can redistribute it and/or modify -! it under the terms of the GNU General Library Public License as published -! by the Free Software Foundation; either version 2 of the License,or -! (at your option) any later version. -! -! MESA 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 Library General Public License for more details. -! -! You should have received a copy of the GNU Library General Public License -! along with this software; if not,write to the Free Software -! Foundation,Inc.,59 Temple Place,Suite 330,Boston,MA 02111-1307 USA -! -! *********************************************************************** - -module test_square_quad - - use mtx_lib - use mtx_def - use utils_lib, only: mesa_error - - implicit none - -contains - - subroutine do_test_square_quad - call test_square_quad1 - call test_square_quad2 - call test_square_quad_inv - end subroutine do_test_square_quad - - subroutine test_square_quad_inv - - integer, parameter :: n = 3, nrhs = 1 - integer :: info, ipiv(n) - real(16) :: A1(n, n), B1(n), A2(n, n), B2(n) - real(16) :: A1_init(n, n), A2_init(n, n), X(n) - - include 'formats' - - write (*, *) 'test_square_quad_inv' - write (*, *) - - A1(1, 1:n) = (/3.14_16, 7.5_16, 0.00_16/) - A1(2, 1:n) = (/4.1_16, 3.2_16, 0.3_16/) - A1(3, 1:n) = (/0.00_16, 1.0_16, 4.1_16/) - A1_init = A1 - - A2(1, 1:n) = (/0.0_16, 3.1_16, 0.0_16/) - A2(2, 1:n) = (/4.7_16, 6.2_16, 0.0_16/) - A2(3, 1:n) = (/3.2_16, 0.0_16, 0.31_16/) - A2_init = A2 - - B1(1:n) = (/1.0_16, 2.0_16, 3.0_16/) - B2(1:n) = (/1.1_16, 2.1_16, 3.1_16/) - - info = 0 - call QGETRF(n, n, A1, n, ipiv, info) - if (info /= 0) then - write (*, *) 'singular matrix?', info - call mesa_error(__FILE__, __LINE__) - end if - X = B1 - ! solve A1*X = B1 - call QGETRS('N', n, nrhs, A1, n, ipiv, X, n, info) - if (info /= 0) call mesa_error(__FILE__, __LINE__) - - write (*, 1) 'B1', B1(1:n) - ! prod = A1_init*X; should get prod == B1 - - write (*, *) 'ipiv1', ipiv(1:n) - write (*, *) - - call QGETRF(n, n, A2, n, ipiv, info) - if (info /= 0) then - write (*, *) 'singular matrix?', info - call mesa_error(__FILE__, __LINE__) - end if - X = B2 - ! solve A2*X = B2 - call QGETRS('N', n, nrhs, A2, n, ipiv, X, n, info) - if (info /= 0) call mesa_error(__FILE__, __LINE__) - - write (*, 1) 'B2', B2(1:n) - - end subroutine test_square_quad_inv - - subroutine test_square_quad2 - - integer, parameter :: n = 3, nrhs = 1 - integer :: info, ipiv(n) - real(16) :: A1(n, n), B1(n, nrhs), A2(n, n), B2(n, nrhs) - - include 'formats' - - write (*, *) 'test test_square_quad2' - write (*, *) - - A1(1, 1:n) = (/3.14_16, 7.5_16, 0.0_16/) - A1(2, 1:n) = (/4.1_16, 3.2_16, 0.3_16/) - A1(3, 1:n) = (/0.00_16, 1.0_16, 4.1_16/) - - A2(1, 1:n) = (/4.7_16, 6.2_16, 0.0_16/) - A2(2, 1:n) = (/3.2_16, 0.0_16, 0.31_16/) - A2(3, 1:n) = (/0.0_16, 3.1_16, 0.0_16/) - - B1(1:n, 1) = (/1.0_16, 2.0_16, 3.0_16/) - B2(1:n, 1) = (/1.1_16, 2.1_16, 3.1_16/) - - info = 0 - call QGETRF(n, n, A1, n, ipiv, info) - if (info /= 0) then - write (*, *) 'singular matrix?', info - call mesa_error(__FILE__, __LINE__) - end if - call QGETRS('N', n, nrhs, A1, n, ipiv, B1, n, info) - if (info /= 0) call mesa_error(__FILE__, __LINE__) - - write (*, 1) 'B1', B1(1:n, 1) - - call QGETRF(n, n, A2, n, ipiv, info) - if (info /= 0) then - write (*, *) 'singular matrix' - call mesa_error(__FILE__, __LINE__) - end if - call QGETRS('N', n, nrhs, A2, n, ipiv, B2, n, info) - if (info /= 0) call mesa_error(__FILE__, __LINE__) - - write (*, 1) 'B2', B2(1:n, 1) - write (*, *) - - end subroutine test_square_quad2 - - subroutine test_square_quad1 - - integer, parameter :: n = 4, nrhs = 1 - integer :: info, ipiv(n) - real(16) :: A(n, n), B(n, nrhs), A2(n, n) - - include 'formats' - - A(1, 1:n) = (/1.80_16, 2.88_16, 2.05_16, 0.00_16/) - A(2, 1:n) = (/5.25_16, -2.95_16, -0.95_16, -3.80_16/) - A(3, 1:n) = (/0.00_16, 0.00_16, -2.90_16, -1.04_16/) - A(4, 1:n) = (/-1.11_16, 0.00_16, -0.59_16, 0.80_16/) - B(1:n, 1) = (/4.35_16, 5.05_16, 3.04_16, -2.05_16/) - - A2 = A - - write (*, *) ' test_square_quad1' - write (*, *) - - info = 0 - call QGETRF(n, n, A, n, ipiv, info) - if (info /= 0) then - write (*, *) 'singular matrix?', info - call mesa_error(__FILE__, __LINE__) - end if - - call QGETRS('N', n, nrhs, A, n, ipiv, B, n, info) - if (info /= 0) call mesa_error(__FILE__, __LINE__) - - write (*, 1) 'B', B(1:n, 1) - write (*, *) - - end subroutine test_square_quad1 - -end module test_square_quad diff --git a/mtx/test/test_output b/mtx/test/test_output index 3e9e193ca..686d99a97 100644 --- a/mtx/test/test_output +++ b/mtx/test/test_output @@ -13,21 +13,6 @@ A1_init*X 1.0000000000000000D+00 2.0000000000000000D+00 2.9999999999999996D+00 B2 1.1000000000000001D+00 2.1000000000000001D+00 3.1000000000000001D+00 A2_init*X 1.1000000000000001D+00 2.1000000000000001D+00 3.1000000000000001D+00 - test_square_quad1 - - B -1.6215408275458856D-01 1.6391349023885566D+00 -3.8454229229650589D-02 -2.8158487838788589D+00 - - test test_square_quad2 - - B1 4.8857961179302775D-01 -7.1218664137347617D-02 7.4907772296032869D-01 - B2 -1.0851063829787234D+00 1.0000000000000000D+00 1.7975291695264242D+01 - - test_square_quad_inv - - B1 1.0000000000000000D+00 2.0000000000000000D+00 3.0000000000000000D+00 - ipiv1 2 2 3 - - B2 1.1000000000000000D+00 2.1000000000000000D+00 3.1000000000000000D+00 testing with nvar,nz 25 20 bcyclic_dble done bcyclic_dble From 73096b3957b8d38256045060a72a679a79a75c57 Mon Sep 17 00:00:00 2001 From: Vincent Vanlaer Date: Thu, 8 Aug 2024 22:19:30 -0400 Subject: [PATCH 2/3] mtx: eliminate all quad precision code None of this code is used in the other parts of MESA. --- mtx/private/mtx_support.f90 | 296 ------------- mtx/public/mtx_decsolblk_quad.dek | 11 - mtx/public/mtx_def.f | 1 - mtx/public/mtx_lib.f90 | 34 +- mtx/public/mtx_null_decsol.dek | 12 - mtx/public/mtx_solve_routines_quad.inc | 544 ------------------------ mtx/test/make/makefile_base | 20 +- mtx/test/src/test_block_tridiagonal.f90 | 137 +----- mtx/test/src/test_mtx.f90 | 2 - 9 files changed, 17 insertions(+), 1040 deletions(-) delete mode 100644 mtx/public/mtx_decsolblk_quad.dek delete mode 100644 mtx/public/mtx_solve_routines_quad.inc diff --git a/mtx/private/mtx_support.f90 b/mtx/private/mtx_support.f90 index 3edf6f51c..53c4ed4b7 100644 --- a/mtx/private/mtx_support.f90 +++ b/mtx/private/mtx_support.f90 @@ -558,39 +558,6 @@ subroutine do_dense_to_col_sparse_0_based_qp( & end subroutine do_dense_to_col_sparse_0_based_qp - subroutine do_quad_dense_to_col_sparse_0_based( & - n,ndim,a,nzmax,nz,colptr,rowind,values,diags,ierr) - integer, intent(in) :: n,ndim,nzmax - real(qp), intent(in) :: a(:,:) ! (ndim,n) - integer, intent(inout) :: colptr(:) ! (n+1) - integer, intent(inout) :: rowind(:) ! (nzmax) - real(qp), intent(out) :: values(:) ! (nzmax) - logical, intent(in) :: diags - integer, intent(out) :: nz,ierr - integer :: i,j - ierr = 0 - nz = 0 - do j=1,n - colptr(j) = nz ! index in values of first entry in this column - do i=1,n - if (a(i,j) == 0) then - if (i /= j) cycle ! not a diagonal - if (.not. diags) cycle - ! else include diagonals even if are == 0 - end if - nz = nz+1 - if (nz > nzmax) then - ierr = j - return - end if - values(nz) = a(i,j) - rowind(nz) = i-1 - end do - end do - colptr(n+1) = nz - end subroutine do_quad_dense_to_col_sparse_0_based - - subroutine do_column_sparse_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr) integer, intent(in) :: n,ndim,nz real(dp), intent(inout) :: a(ndim,n) @@ -617,28 +584,6 @@ subroutine do_column_sparse_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr) end subroutine do_column_sparse_to_dense - subroutine do_quad_column_sparse_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr) - integer, intent(in) :: n,ndim,nz - real(qp), intent(out) :: a(ndim,n) - integer, intent(in) :: colptr(n+1),rowind(nz) - real(qp), intent(in) :: values(nz) - integer, intent(out) :: ierr - integer :: i,j,k - ierr = 0 - a = 0 - do j=1,n - do k=colptr(j),colptr(j+1)-1 - i = rowind(k) - if (i > n) then - ierr = j - return - endif - a(i,j) = values(k) - end do - end do - end subroutine do_quad_column_sparse_to_dense - - subroutine do_col_sparse_0_based_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr) integer, intent(in) :: n,ndim,nz real(dp), intent(inout) :: a(ndim,n) @@ -685,56 +630,6 @@ subroutine do_block_dble_mv(nvar, nz, lblk, dblk, ublk, b, prod) end subroutine do_block_dble_mv - subroutine do_block_mv_quad(lblk, dblk, ublk, b, prod) - real(qp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz) - real(qp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz) - real(qp), pointer, dimension(:,:), intent(inout) :: prod ! (nvar,nz) - integer :: nvar, nz, k - include 'formats' - nvar = size(b,dim=1) - nz = size(b,dim=2) -!$OMP PARALLEL DO PRIVATE(k) - do k = 1, nz - prod(:,k) = 0 - call qdgemv(nvar,nvar,dblk(:,:,k),nvar,b(:,k),prod(:,k)) - if (k > 1) then - call qdgemv(nvar,nvar,lblk(:,:,k),nvar,b(:,k-1),prod(:,k)) - end if - if (k < nz) then - call qdgemv(nvar,nvar,ublk(:,:,k),nvar,b(:,k+1),prod(:,k)) - end if - end do -!$OMP END PARALLEL DO - - contains - - subroutine qdgemv(m,n,a,lda,x,y) - ! y := alpha*a*x + beta*y - use const_def, only: dp - - integer lda,m,n - real(qp) a(lda,*),x(*),y(*) - real(qp) :: tmp - ! trans = 'n' - ! alpha = 1 - ! beta = 1 - ! incx = 1 - ! incy = 1 - integer :: j, i - do j = 1,n - tmp = x(j) - if (tmp.ne.0d0) then - do i = 1,m - y(i) = y(i) + tmp*a(i,j) - end do - end if - end do - end subroutine qdgemv - - - end subroutine do_block_mv_quad - - subroutine do_LU_factored_block_dble_mv(lblk, dblk, ublk, b, ipiv, prod) real(dp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz) real(dp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz) @@ -836,21 +731,6 @@ subroutine do_multiply_xa(n, A1, x, b) end subroutine do_multiply_xa - subroutine do_quad_multiply_xa(n, A1, x, b) - ! calculates b = x*A - integer, intent(in) :: n - real(qp), pointer, intent(in) :: A1(:) ! =(n, n) - real(qp), pointer, intent(in) :: x(:) ! (n) - real(qp), pointer, intent(inout) :: b(:) ! (n) - integer :: j - real(qp), pointer :: A(:,:) ! (n, n) - A(1:n,1:n) => A1(1:n*n) - do j = 1, n - b(j) = dot_product(x(1:n),A(1:n,j)) - end do - end subroutine do_quad_multiply_xa - - subroutine do_multiply_xa_plus_c(n, A1, x, c, b) ! calculates b = x*A + c integer, intent(in) :: n @@ -867,22 +747,6 @@ subroutine do_multiply_xa_plus_c(n, A1, x, c, b) end subroutine do_multiply_xa_plus_c - subroutine do_quad_multiply_xa_plus_c(n, A1, x, c, b) - ! calculates b = x*A + c - integer, intent(in) :: n - real(qp), pointer, intent(in) :: A1(:) ! =(n, n) - real(qp), pointer, intent(in) :: x(:) ! (n) - real(qp), pointer, intent(in) :: c(:) ! (n) - real(qp), pointer, intent(inout) :: b(:) ! (n) - integer :: j - real(qp), pointer :: A(:,:) ! (n, n) - A(1:n,1:n) => A1(1:n*n) - do j = 1, n - b(j) = dot_product(x(1:n),A(1:n,j)) + c(j) - end do - end subroutine do_quad_multiply_xa_plus_c - - subroutine do_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1) ! calculates b = x*A integer, intent(in) :: nvar, nz @@ -919,42 +783,6 @@ subroutine do_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1) end subroutine do_block_multiply_xa - subroutine do_quad_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1) - ! calculates b = x*A - integer, intent(in) :: nvar, nz - real(qp), dimension(:), intent(in), pointer :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nz) - real(qp), intent(in), pointer :: x1(:) ! =(nvar,nz) - real(qp), intent(inout), pointer :: b1(:) ! =(nvar,nz) - integer :: k, shift, shift2, nvar2 - real(qp), pointer, dimension(:) :: p1, p2, p3, p4 - nvar2 = nvar*nvar -!$OMP PARALLEL DO PRIVATE(k,shift,shift2,p1,p2,p3,p4) - do k = 1, nz - shift = nvar*(k-1) - shift2 = nvar2*(k-1) - p1(1:nvar2) => dblk1(shift2+1:shift2+nvar2) - p2(1:nvar) => x1(shift+1:shift+nvar) - p3(1:nvar) => b1(shift+1:shift+nvar) - call do_quad_multiply_xa(nvar,p1,p2,p3) - if (k > 1) then - p1(1:nvar2) => ublk1(shift2+1:shift2+nvar2) - p2(1:nvar) => x1(shift+1:shift+nvar) - p3(1:nvar) => b1(shift+1:shift+nvar) - p4(1:nvar) => b1(shift+1:shift+nvar) - call do_quad_multiply_xa_plus_c(nvar,p1,p2,p3,p4) - end if - if (k < nz) then - p1(1:nvar2) => lblk1(shift2+1:shift2+nvar2) - p2(1:nvar) => x1(shift+1+nvar:shift+2*nvar) - p3(1:nvar) => b1(shift+1:shift+nvar) - p4(1:nvar) => b1(shift+1:shift+nvar) - call do_quad_multiply_xa_plus_c(nvar,p1,p2,p3,p4) - end if - end do -!$OMP END PARALLEL DO - end subroutine do_quad_block_multiply_xa - - subroutine do_band_multiply_xa(n, kl, ku, ab1, ldab, x, b) ! calculates b = x*a = transpose(a)*x integer, intent(in) :: n @@ -988,38 +816,6 @@ subroutine do_band_multiply_xa(n, kl, ku, ab1, ldab, x, b) end subroutine do_band_multiply_xa - subroutine do_quad_band_multiply_xa(n, kl, ku, ab, ldab, x, b) - ! calculates b = x*a = transpose(a)*x - integer, intent(in) :: n - ! the number of linear equations, i.e., the order of the - ! matrix a. n >= 0. - integer, intent(in) :: kl - ! the number of subdiagonals within the band of a. kl >= 0. - integer, intent(in) :: ku - ! the number of superdiagonals within the band of a. ku >= 0. - integer, intent(in) :: ldab - ! the leading dimension of the array ab. ldab >= kl+ku+1. - real(qp), intent(in) :: ab(:,:) ! (ldab, n) - ! the matrix a in band storage, in rows 1 to kl+ku+1; - ! the j-th column of a is stored in the j-th column of the - ! array ab as follows: - ! ab(ku+1+i-j, j) = a(i, j) for max(1, j-ku)<=i<=min(n, j+kl) - real(qp), intent(in) :: x(:) ! (n) - ! the input vector to be multiplied by the matrix. - real(qp), intent(inout) :: b(:) ! (n) - ! on exit, set to matrix product of x*a = b - integer :: i, j, k - do j = 1, n - k = ku+1-j - b(j) = 0 - do i = max(1, j-ku), min(n, j+kl) - b(j) = b(j) + x(i)*ab(k+i, j) - end do - end do - end subroutine do_quad_band_multiply_xa - - - subroutine do_clip_blocks( & mblk, clip_limit, lmat, dmat, umat, dmat_nnz, total_nnz) integer, intent(in) :: mblk @@ -1106,52 +902,6 @@ subroutine read_hbcode1(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr end subroutine read_hbcode1 - subroutine read_hbcode1_quad(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr) - - CHARACTER TITLE*72 , KEY*8 , MXTYPE*3 , & - PTRFMT*16, INDFMT*16, VALFMT*20, RHSFMT*20 - - INTEGER TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, & - iounit, NROW , NCOL , NNZERO, NELTVL - - INTEGER, pointer :: COLPTR (:), ROWIND (:) - - REAL(qp), pointer :: VALUES (:) - integer, intent(out) :: ierr - - integer :: i - ierr = 0 - READ (iounit, 1000, iostat=ierr ) TITLE , KEY , & - TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, & - MXTYPE, NROW , NCOL , NNZERO, NELTVL, & - PTRFMT, INDFMT, VALFMT, RHSFMT - if (ierr /= 0) return - 1000 FORMAT ( A72, A8 / 5I14 / A3, 11X, 4I14 / 2A16, 2A20 ) - - allocate(VALUES(NNZERO), ROWIND(NNZERO), COLPTR(NCOL+1), stat=ierr) - if (ierr /= 0) return - - READ (iounit, PTRFMT, iostat=ierr ) ( COLPTR (I), I = 1, NCOL+1 ) - if (ierr /= 0) return - - READ (iounit, INDFMT, iostat=ierr ) ( ROWIND (I), I = 1, NNZERO ) - if (ierr /= 0) return - - IF ( VALCRD .GT. 0 ) THEN - -! ---------------------- -! ... READ MATRIX VALUES -! ---------------------- - - READ (iounit, VALFMT, iostat=ierr ) ( VALUES (I), I = 1, NNZERO ) - if (ierr /= 0) return - - ENDIF - - - end subroutine read_hbcode1_quad - - subroutine write_hbcode1(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr) CHARACTER TITLE*72 , KEY*8 , MXTYPE*3 , & @@ -1304,37 +1054,6 @@ subroutine read_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr) end subroutine read_block_tridiagonal - subroutine read_quad_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr) - integer, intent(in) :: iounit - integer, intent(out) :: nvar, nblk - real(qp), pointer, dimension(:) :: lblk1,dblk1,ublk1 ! will be allocated - integer, intent(out) :: ierr - integer :: k - real(qp), pointer, dimension(:,:,:) :: lblk,dblk,ublk - ierr = 0 - read(iounit,*,iostat=ierr) nvar, nblk - if (ierr /= 0) return - allocate(lblk1(nvar*nvar*nblk), dblk1(nvar*nvar*nblk), ublk1(nvar*nvar*nblk), stat=ierr) - if (ierr /= 0) return - lblk(1:nvar,1:nvar,1:nblk) => lblk1(1:nvar*nvar*nblk) - dblk(1:nvar,1:nvar,1:nblk) => dblk1(1:nvar*nvar*nblk) - ublk(1:nvar,1:nvar,1:nblk) => ublk1(1:nvar*nvar*nblk) - do k=1,nblk - if (k > 1) then - call read1_sparse_block_quad(iounit, nvar, lblk(:,:,k), ierr) - if (ierr /= 0) return - end if - call read1_sparse_block_quad(iounit, nvar, dblk(:,:,k), ierr) - if (ierr /= 0) return - if (k < nblk) then - call read1_sparse_block_quad(iounit, nvar, ublk(:,:,k), ierr) - if (ierr /= 0) return - end if - end do - - end subroutine read_quad_block_tridiagonal - - subroutine read1_sparse_block(iounit, nvar, blk, ierr) integer, intent(in) :: iounit, nvar real(dp) :: blk(:,:) ! (nvar,nvar) @@ -1350,21 +1069,6 @@ subroutine read1_sparse_block(iounit, nvar, blk, ierr) end subroutine read1_sparse_block - subroutine read1_sparse_block_quad(iounit, nvar, blk, ierr) - integer, intent(in) :: iounit, nvar - real(qp) :: blk(:,:) ! (nvar,nvar) - integer, intent(out) :: ierr - integer :: nnz, nrow, ncol - integer, pointer :: rowind(:), colptr(:) - real(qp), pointer :: values(:) - ierr = 0 - call read_hbcode1_quad(iounit, nrow, ncol, nnz, values, rowind, colptr,ierr) - if (ierr /= 0 .or. nrow /= nvar .or. nrow /= ncol) return - call do_quad_column_sparse_to_dense(nrow,ncol,blk,nnz,colptr,rowind,values,ierr) - deallocate(colptr,rowind,values) - end subroutine read1_sparse_block_quad - - subroutine write_block_tridiagonal(iounit,nvar,nblk,lblk,dblk,ublk,ierr) integer, intent(in) :: iounit, nvar, nblk real(dp), intent(in), dimension(:,:,:) :: lblk,dblk,ublk diff --git a/mtx/public/mtx_decsolblk_quad.dek b/mtx/public/mtx_decsolblk_quad.dek deleted file mode 100644 index 3d40a9355..000000000 --- a/mtx/public/mtx_decsolblk_quad.dek +++ /dev/null @@ -1,11 +0,0 @@ - ! block tridiagonal - subroutine decsolblk_quad(iop,caller_id,nvar,nz,lblk,dblk,ublk,brhs,ipiv,lrd,rpar_decsol,lid,ipar_decsol,ierr) - use const_def, only: qp, dp - integer, intent(in) :: caller_id, nvar, nz, iop, lrd, lid - integer, pointer, intent(inout) :: ipiv(:) ! =(nvar,nz) - real(qp), dimension(:), pointer, intent(inout) :: lblk, dblk, ublk ! =(nvar,nvar,nz) - real(qp), pointer, intent(inout) :: brhs(:) ! =(nvar,nz) - real(dp), pointer, intent(inout) :: rpar_decsol(:) ! (lrd) - integer, pointer, intent(inout) :: ipar_decsol(:) ! (lid) - integer, intent(out) :: ierr - end subroutine decsolblk_quad diff --git a/mtx/public/mtx_def.f b/mtx/public/mtx_def.f index a00ba566c..6e0858849 100644 --- a/mtx/public/mtx_def.f +++ b/mtx/public/mtx_def.f @@ -30,7 +30,6 @@ module mtx_def ! matrix solver options integer, parameter :: lapack = 1 integer, parameter :: block_thomas_dble = 2 - integer, parameter :: block_thomas_quad = 3 integer, parameter :: block_thomas_refine = 4 integer, parameter :: bcyclic_dble = 5 diff --git a/mtx/public/mtx_lib.f90 b/mtx/public/mtx_lib.f90 index d72b8c5fb..002314f80 100644 --- a/mtx/public/mtx_lib.f90 +++ b/mtx/public/mtx_lib.f90 @@ -26,7 +26,7 @@ module mtx_lib - use const_def, only: dp, qp + use const_def, only: dp implicit none @@ -237,15 +237,6 @@ subroutine mtx_read_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr) call read_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr) end subroutine mtx_read_block_tridiagonal - subroutine mtx_read_quad_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr) - use mtx_support, only: read_quad_block_tridiagonal - integer, intent(in) :: iounit - integer, intent(out) :: nvar, nblk - real(qp), pointer, dimension(:) :: lblk1,dblk1,ublk1 ! =(nvar,nvar,nblk) will be allocated - integer, intent(out) :: ierr - call read_quad_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr) - end subroutine mtx_read_quad_block_tridiagonal - ! BCYCLIC multi-thread block tridiagonal include "mtx_bcyclic_dble_decsol.dek" ! S.P.Hirshman, K.S.Perumalla, V.E.Lynch, & R.Sanchez, @@ -264,17 +255,6 @@ subroutine block_dble_mv(nvar, nz, lblk, dblk, ublk, b, prod) end subroutine block_dble_mv - subroutine block_quad_mv(lblk, dblk, ublk, b, prod) - ! set prod = A*b with A = block tridiagonal given by lblk, dblk, ublk - use mtx_support, only: do_block_mv_quad - real(qp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz) - real(qp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz) - real(qp), pointer, dimension(:,:), intent(inout) :: prod ! (nvar,nz) - call do_block_mv_quad(lblk, dblk, ublk, b, prod) - end subroutine block_quad_mv - - - subroutine multiply_xa(n, A1, x, b) ! calculates b = x*A use mtx_support, only: do_multiply_xa @@ -296,18 +276,6 @@ subroutine block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1) call do_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1) end subroutine block_multiply_xa - - subroutine quad_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1) - ! calculates b = x*A - use mtx_support, only: do_quad_block_multiply_xa - integer, intent(in) :: nvar, nz - real(qp), dimension(:), intent(in), pointer :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nz) - real(qp), intent(in), pointer :: x1(:) ! =(nvar,nz) - real(qp), intent(inout), pointer :: b1(:) ! =(nvar,nz) - call do_quad_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1) - end subroutine quad_block_multiply_xa - - subroutine band_multiply_xa(n, kl, ku, ab1, ldab, x, b) ! calculates b = x*a = transpose(a)*x use mtx_support, only: do_band_multiply_xa diff --git a/mtx/public/mtx_null_decsol.dek b/mtx/public/mtx_null_decsol.dek index 103989c26..4b45f1788 100644 --- a/mtx/public/mtx_null_decsol.dek +++ b/mtx/public/mtx_null_decsol.dek @@ -22,18 +22,6 @@ end subroutine null_decsolblk - subroutine null_decsolblk_quad(iop,caller_id,nvar,nz,lblk,dblk,ublk,brhs,ipiv,lrd,rpar_decsol,lid,ipar_decsol,ierr) - integer, intent(in) :: caller_id, iop, nvar, nz, lrd, lid - integer, pointer, intent(inout) :: ipiv(:) ! =(nvar,nz) - real(qp), dimension(:), pointer, intent(inout) :: lblk, dblk, ublk ! =(nvar,nvar,nz) - real(qp), pointer, intent(inout) :: brhs(:) ! =(nvar,nz) - real(dp), pointer, intent(inout) :: rpar_decsol(:) ! (lrd) - integer, pointer, intent(inout) :: ipar_decsol(:) ! (lid) - integer, intent(out) :: ierr - ierr = -1 - end subroutine null_decsolblk_quad - - subroutine null_decsolc(iop,n,ndim,ar,ai,ml,mu,br,bi,ip,lcd,cpar_decsol,lrd,rpar_decsol,lid,ipar_decsol,ier) integer, intent(in) :: iop, n, ndim, lcd, lrd, lid, ml, mu integer, intent(inout) :: ip(n) diff --git a/mtx/public/mtx_solve_routines_quad.inc b/mtx/public/mtx_solve_routines_quad.inc deleted file mode 100644 index 983af81c8..000000000 --- a/mtx/public/mtx_solve_routines_quad.inc +++ /dev/null @@ -1,544 +0,0 @@ - - - subroutine my_getf2_qp(m, a, lda, ipiv, info) - integer :: info, lda, m - integer :: ipiv(:) - real(qp) :: a(:,:) - real(qp), parameter :: one=1, zero=0 - integer :: i, j, jp, ii, jj, n, mm - real(qp) :: tmp, da - do j = 1, m - info = 0 - jp = j - 1 + maxloc(abs(a(j:lda,j)),dim=1) - ipiv( j ) = jp - if( a( jp, j ).ne.zero ) then - if( jp.ne.j ) then ! swap a(j,:) and a(jp,:) - !!!! not with quad !xxdec$ simd private(tmp) - do i=1,m - tmp = a(j,i) - a(j,i) = a(jp,i) - a(jp,i) = tmp - end do - end if - if( j.lt.m ) then - if (a(j,j) == 0) then - info = -1 - return - end if - da = one / a( j, j ) - ! not with quad !xxdec$ simd - do i = 1, m-j - a( j+i, j ) = da*a( j+i, j ) - end do - end if - else if( info.eq.0 ) then - info = j - end if - if( j.lt.m ) then - !call dger( m-j, m-j, -one, a( j+1, j ), 1, a( j, j+1 ), lda, a( j+1, j+1 ), lda ) - do jj = j+1, m - ! not with quad !xxdec$ simd - do ii = j+1, m - a(ii,jj) = a(ii,jj) - a(ii,j)*a(j,jj) - end do - end do - end if - end do - end subroutine my_getf2_qp - - - subroutine my_getrs_qp( n, nrhs, a, lda, ipiv, b, ldb, info ) - integer :: info, lda, ldb, n, nrhs - integer :: ipiv(:) - real(qp) :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * ) - real(qp), parameter :: one=1, zero=0 - real(qp) :: temp - integer :: i,j,k, n32, ix, ip - info = 0 - call my_laswp_qp(nrhs, b, ldb, 1, n, ipiv, 1 ) - !call dtrsm( 'left', 'lower', 'no transpose', 'unit', n, nrhs, one, a, lda, b, ldb ) - do j = 1,nrhs - do k = 1,n - if (b(k,j).ne.zero) then - ! not with quad !xxdec$ simd - do i = k + 1,n - b(i,j) = b(i,j) - b(k,j)*a(i,k) - end do - end if - end do - end do - !call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, nrhs, one, a, lda, b, ldb ) - do j = 1,nrhs - do k = n,1,-1 - if (b(k,j).ne.zero) then - if (a(k,k) == 0) then - info = -1 - return - end if - b(k,j) = b(k,j)/a(k,k) - ! not with quad !xxdec$ simd - do i = 1,k - 1 - b(i,j) = b(i,j) - b(k,j)*a(i,k) - end do - end if - end do - end do - - end subroutine my_getrs_qp - - - subroutine my_laswp_qp( n, a, lda, k1, k2, ipiv, incx ) - integer :: incx, k1, k2, lda, n - integer :: ipiv(:) - real(qp) :: a(:,:) ! a( lda, * ) - integer :: i, i1, i2, inc, ip, ix, ix0, j, k, n32 - real(qp) :: temp - ! interchange row i with row ipiv(i) for each of rows k1 through k2. - if( incx.gt.0 ) then - ix0 = k1 - i1 = k1 - i2 = k2 - inc = 1 - else if( incx.lt.0 ) then - ix0 = 1 + ( 1-k2 )*incx - i1 = k2 - i2 = k1 - inc = -1 - else - return - end if - ix = ix0 - do i = i1, i2, inc - ip = ipiv( ix ) - if( ip.ne.i ) then - !!!! not with quad !xxdec$ simd private(temp) - do k = 1, n - temp = a( i, k ) - a( i, k ) = a( ip, k ) - a( ip, k ) = temp - end do - end if - ix = ix + incx - end do - end subroutine my_laswp_qp - - - subroutine my_getrs1_qp( n, a, lda, ipiv, b, ldb, info ) - integer :: info, lda, ldb, n - integer :: ipiv(:) - real(qp) :: a(:,:), b(:) - real(qp), parameter :: one=1, zero=0 - real(qp) :: temp - integer :: i,k - info = 0 - !!!! not with quad !xxdec$ simd private(temp) ! ?? <<<< buggy in ifort 14.0.3 - do i = 1, n - temp = b(i) - b(i) = b(ipiv(i)) - b(ipiv(i)) = temp - end do - do k = 1,n - if (b(k).ne.zero) then - ! not with quad !xxdec$ simd - do i = k + 1,n - b(i) = b(i) - b(k)*a(i,k) - end do - end if - end do - do k = n,1,-1 - if (b(k).ne.zero) then - b(k) = b(k)/a(k,k) - ! not with quad !xxdec$ simd - do i = 1,k - 1 - b(i) = b(i) - b(k)*a(i,k) - end do - end if - end do - end subroutine my_getrs1_qp - - - subroutine my_gemm0_p1_qp(m,n,k,a,lda,b,ldb,c,ldc) ! c := -a*b - integer, intent(in) :: k,lda,ldb,ldc,m,n - real(qp), dimension(:,:) :: a, b, c ! a(lda,*),b(ldb,*),c(ldc,*) - integer :: j, i - real(qp), parameter :: zero=0 - include 'formats' - ! transa = 'n' - ! transb = 'n' - ! alpha = -1 - ! beta = 0 - ! assumes other args are valid - c(1:m,1:n) = zero - call my_gemm_p1_qp(m,n,k,a,lda,b,ldb,c,ldc) - end subroutine my_gemm0_p1_qp - - - subroutine my_gemm_plus_mm_qp(m,n,k,a,b,d,e,c) ! c := c + a*b + d*e - integer, intent(in) :: k,m,n - real(qp), dimension(:,:) :: a, b, c, d, e - real(qp) :: tmp_b, tmp_e - real(qp), parameter :: zero=0 - integer :: j, i, l - do j = 1,n - do l = 1,k - tmp_b = b(l,j) - tmp_e = e(l,j) - if (tmp_b .ne. zero) then - if (tmp_e .ne. zero) then - ! not with quad !xxdec$ simd - do i = 1,m - c(i,j) = c(i,j) + tmp_b*a(i,l) + tmp_e*d(i,l) - end do - else - ! not with quad !xxdec$ simd - do i = 1,m - c(i,j) = c(i,j) + tmp_b*a(i,l) - end do - end if - else if (tmp_e .ne. zero) then - ! not with quad !xxdec$ simd - do i = 1,m - c(i,j) = c(i,j) + tmp_e*d(i,l) - end do - end if - end do - end do - end subroutine my_gemm_plus_mm_qp - - - subroutine my_gemm_qp(m,n,k,a,lda,b,ldb,c,ldc) ! c := c - a*b - integer, intent(in) :: k,lda,ldb,ldc,m,n - real(qp), dimension(:,:) :: a, b, c ! a(lda,*),b(ldb,*),c(ldc,*) - real(qp) :: tmp - real(qp), parameter :: zero=0 - integer :: j, i, l - ! transa = 'n' - ! transb = 'n' - ! alpha = -1 - ! beta = 1 - ! assumes other args are valid - do j = 1,n - do l = 1,k - tmp = b(l,j) - if (tmp .ne. zero) then - ! not with quad !xxdec$ simd - do i = 1,m - c(i,j) = c(i,j) - tmp*a(i,l) - end do - end if - end do - end do - end subroutine my_gemm_qp - - - subroutine my_gemm_p1_qp(m,n,k,a,lda,b,ldb,c,ldc) ! c := c + a*b - integer, intent(in) :: k,lda,ldb,ldc,m,n - real(qp), dimension(:,:) :: a, b, c ! a(lda,*),b(ldb,*),c(ldc,*) - real(qp) :: tmp - real(qp), parameter :: zero=0 - integer :: j, i, l - ! transa = 'n' - ! transb = 'n' - ! alpha = 1 - ! beta = 1 - ! assumes other args are valid - do j = 1,n - do l = 1,k - tmp = b(l,j) - if (tmp .ne. zero) then - ! not with quad !xxdec$ simd - do i = 1,m - c(i,j) = c(i,j) + tmp*a(i,l) - end do - end if - end do - end do - end subroutine my_gemm_p1_qp - - - subroutine my_gemv_mv_qp(m,n,a,x,b,z,y) ! y = y - a*x - b*z - integer lda,m,n - real(qp) :: a(:,:), b(:,:) - real(qp) :: x(:), z(:), y(:) - real(qp) :: tmp_x, tmp_z - real(qp), parameter :: zero=0 - integer :: j, i - do j = 1,n - tmp_x = x(j) - tmp_z = z(j) - if (tmp_x.ne.zero) then - if (tmp_z /= zero) then - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) - tmp_x*a(i,j) - tmp_z*b(i,j) - end do - else - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) - tmp_x*a(i,j) - end do - end if - else if (tmp_z /= zero) then - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) - tmp_z*b(i,j) - end do - end if - end do - end subroutine my_gemv_mv_qp - - - subroutine my_gemv_qp(m,n,a,lda,x,y) ! y = y - a*x - integer lda,m,n - real(qp) :: a(:,:) ! (lda,*) - real(qp) :: x(:), y(:) - real(qp) :: tmp - real(qp), parameter :: zero=0 - ! trans = 'n' - ! alpha = -1 - ! beta = 1 - ! incx = 1 - ! incy = 1 - integer :: j, i - do j = 1,n - tmp = x(j) - if (tmp.ne.zero) then - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) - tmp*a(i,j) - end do - end if - end do - end subroutine my_gemv_qp - - - subroutine my_gemv_p_mv_qp(m,n,a,x,b,z,y) ! y = y + a*x + b*z - integer lda,m,n - real(qp) :: a(:,:), b(:,:) - real(qp) :: x(:), z(:), y(:) - real(qp) :: tmp_x, tmp_z - real(qp), parameter :: zero=0 - integer :: j, i - do j = 1,n - tmp_x = x(j) - tmp_z = z(j) - if (tmp_x.ne.zero) then - if (tmp_z /= zero) then - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) + tmp_x*a(i,j) + tmp_z*b(i,j) - end do - else - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) + tmp_x*a(i,j) - end do - end if - else if (tmp_z /= zero) then - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) + tmp_z*b(i,j) - end do - end if - end do - end subroutine my_gemv_p_mv_qp - - - subroutine my_gemv_p1_qp(m,n,a,lda,x,y) ! y = y + a*x - integer lda,m,n - real(qp) :: a(:,:) ! (lda,*) - real(qp) :: x(:), y(:) - real(qp) :: tmp - real(qp), parameter :: zero=0 - ! trans = 'n' - ! alpha = -1 - ! beta = 1 - ! incx = 1 - ! incy = 1 - integer :: j, i - do j = 1,n - tmp = x(j) - if (tmp.ne.zero) then - ! not with quad !xxdec$ simd - do i = 1,m - y(i) = y(i) + tmp*a(i,j) - end do - end if - end do - end subroutine my_gemv_p1_qp - - - subroutine my_getf2_no_pivot_qp(m, a, lda, info) - integer :: info, lda, m - real(qp) :: a(:,:) - real(qp), parameter :: one=1, zero=0 - integer :: i, j, ii, jj, n, mm - real(qp) :: tmp, da - do j = 1, m - info = 0 - if( a( j, j ).ne.zero ) then - if( j.lt.m ) then - da = one / a( j, j ) - ! not with quad !xxdec$ simd - do i = 1, m-j - a( j+i, j ) = da*a( j+i, j ) - end do - end if - else if( info.eq.0 ) then - info = j - end if - if( j.lt.m ) then - do jj = j+1, m - ! not with quad !xxdec$ simd - do ii = j+1, m - a(ii,jj) = a(ii,jj) - a(ii,j)*a(j,jj) - end do - end do - end if - end do - end subroutine my_getf2_no_pivot_qp - - - subroutine my_getrs_no_pivot_qp( n, nrhs, a, lda, b, ldb, info ) - integer :: info, lda, ldb, n, nrhs - real(qp) :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * ) - real(qp), parameter :: one=1, zero=0 - real(qp) :: temp - integer :: i,j,k, n32, ix, ip - info = 0 - do j = 1,nrhs - do k = 1,n - if (b(k,j).ne.zero) then - ! not with quad !xxdec$ simd - do i = k + 1,n - b(i,j) = b(i,j) - b(k,j)*a(i,k) - end do - end if - end do - end do - do j = 1,nrhs - do k = n,1,-1 - if (b(k,j).ne.zero) then - b(k,j) = b(k,j)/a(k,k) - ! not with quad !xxdec$ simd - do i = 1,k - 1 - b(i,j) = b(i,j) - b(k,j)*a(i,k) - end do - end if - end do - end do - end subroutine my_getrs_no_pivot_qp - - - subroutine my_getrs1_no_pivot_qp( n, a, lda, b, ldb, info ) - integer :: info, lda, ldb, n - real(qp) :: a(:,:), b(:) - real(qp), parameter :: one=1, zero=0 - real(qp) :: temp - integer :: i,k - info = 0 - do k = 1,n - if (b(k).ne.zero) then - ! not with quad !xxdec$ simd - do i = k + 1,n - b(i) = b(i) - b(k)*a(i,k) - end do - end if - end do - do k = n,1,-1 - if (b(k).ne.zero) then - b(k) = b(k)/a(k,k) - ! not with quad !xxdec$ simd - do i = 1,k - 1 - b(i) = b(i) - b(k)*a(i,k) - end do - end if - end do - end subroutine my_getrs1_no_pivot_qp - - - - subroutine qrdcmp_qp(a,n,c,d,sing) - ! constructs the qr decomposition of a(1:n,1:n). - ! the upper triangular matrix r is returned in the upper triangle of a, except - ! for the diagonal which is returned in d(1:n). the orthogonal matrix q is - ! represented as a product of n-1 householder matrices q_1....q_{n-1}, where - ! q_j = i - u_j * u_j/c_j. the ith component of u_j is zero for i=1..j-1, - ! while the nonzero components are returned in a(i,j) for i=j....n. sing - ! returns as true if singularity is encountered during the decomposition, - ! but the decomposition is still completed. - integer, intent(in) :: n - real(qp), intent(inout) :: a(:,:) - real(qp), intent(out) :: c(:),d(:) - logical, intent(out) :: sing - integer :: i,j,k - real(qp) :: scale,sigma,sum,tau,dum - sing = .false. - do k=1,n-1 - scale = 0.0d0 - do i=k,n - scale = max(scale,abs(a(i,k))) - enddo - if (scale .eq. 0.0) then - sing = .true. - c(k) = 0.0d0 - d(k) = 0.0d0 - else - dum = 1.0d0/scale - do i=k,n - a(i,k) = a(i,k) * dum - enddo - sum = 0.0d0 - do i=k,n - sum = sum + a(i,k)*a(i,k) - enddo - sigma = sign(sqrt(sum),a(k,k)) - a(k,k) = a(k,k) + sigma - c(k) = sigma * a(k,k) - d(k) = -scale * sigma - do j=k+1,n - sum = 0.0d0 - do i=k,n - sum = sum + a(i,k)*a(i,j) - enddo - tau = sum/c(k) - do i=k,n - a(i,j) = a(i,j) - tau * a(i,k) - enddo - enddo - end if - enddo - d(n) = a(n,n) - if (d(n) .eq. 0.0) sing = .true. - end subroutine qrdcmp_qp - - - subroutine qrsolv_qp(a,n,c,d,b) - ! solves the linear system a*x=b, given the qr decomposition from - ! routine qrdcmp. b is input as the right hand side and returns with the - ! solution x. - integer, intent(in) :: n - real(qp), intent(in) :: a(:,:),c(:),d(:) - real(qp), intent(inout) :: b(:) - integer :: i,j - real(qp) :: sum,tau - do j=1,n-1 - sum = 0.0d0 - do i=j,n - sum = sum + a(i,j)*b(i) - enddo - tau = sum/c(j) - do i=j,n - b(i) = b(i) - tau*a(i,j) - enddo - enddo - b(n) = b(n)/d(n) - do i=n-1,1,-1 - sum = 0.0d0 - do j=i+1,n - sum = sum + a(i,j)*b(j) - enddo - b(i) = (b(i) - sum)/d(i) - enddo - end subroutine qrsolv_qp - diff --git a/mtx/test/make/makefile_base b/mtx/test/make/makefile_base index 356e90779..15e335478 100644 --- a/mtx/test/make/makefile_base +++ b/mtx/test/make/makefile_base @@ -9,7 +9,7 @@ MESA_DIR = ../../.. include $(MESA_DIR)/utils/makefile_header SRCS = \ - test_mtx_support.f90 test_block_tri_dble.f90 test_block_tri_quad.f90 \ + test_mtx_support.f90 test_block_tridiagonal.f90 \ test_square.f90 test_mtx.f90 ################################################################# @@ -47,7 +47,7 @@ else endif clean: - -@rm -f *.o *.mod .depend $(TEST) test_block_tri_dble.f90 test_block_tri_quad.f90 + -@rm -f *.o *.mod .depend $(TEST) nodeps : $(.DEFAULT_GOAL) @@ -57,22 +57,6 @@ nodeps : $(.DEFAULT_GOAL) PREPROCESS = $(FC) $(FCfree) $(FC_free_preprocess) -E -test_block_tri_dble.f90: test_block_tridiagonal.f90 -ifneq ($(QUIET),) - @echo PREPROCESS $< - @$(PREPROCESS) -E -DDBLE $< > $@ -else - $(PREPROCESS) -E -DDBLE $< > $@ -endif - -test_block_tri_quad.f90: test_block_tridiagonal.f90 -ifneq ($(QUIET),) - @echo PREPROCESS $< - @$(PREPROCESS) -E $< > $@ -else - $(PREPROCESS) -E $< > $@ -endif - %.o: %.f90 ifneq ($(QUIET),) @echo TEST_COMPILE $< diff --git a/mtx/test/src/test_block_tridiagonal.f90 b/mtx/test/src/test_block_tridiagonal.f90 index 853295b16..3c77bd458 100644 --- a/mtx/test/src/test_block_tridiagonal.f90 +++ b/mtx/test/src/test_block_tridiagonal.f90 @@ -20,59 +20,33 @@ ! ! *********************************************************************** -#ifdef DBLE module test_block_tri_dble -#else - module test_block_tri_quad -#endif - use mtx_lib use mtx_def -#ifdef DBLE use const_def, only: dp -#else - use const_def, only: qp -#endif use utils_lib, only: is_bad, mesa_error implicit none -#ifdef DBLE integer, parameter :: fltp = dp -#else - integer, parameter :: fltp = qp -#endif - integer, parameter :: caller_id = 0 contains -#ifdef DBLE - subroutine do_test_block_tri_dble call test_block(bcyclic_dble, .true.) call test_block(block_thomas_dble, .true.) end subroutine do_test_block_tri_dble -#else - - subroutine do_test_block_tri_quad - call test_block(block_thomas_quad, .true.) - end subroutine do_test_block_tri_quad - -#endif - subroutine test_block(which_decsol_option, for_release) integer, intent(in) :: which_decsol_option logical, intent(in) :: for_release - integer :: nrep real(fltp), pointer :: lblk1(:), dblk1(:), ublk1(:) ! (nvar,nvar,nz) real(fltp), pointer :: lblk(:, :, :), dblk(:, :, :), ublk(:, :, :) ! (nvar,nvar,nz) - real(fltp), pointer :: l_init(:, :, :), d_init(:, :, :), u_init(:, :, :) ! (nvar,nvar,nz) real(fltp), pointer :: x(:, :), xcorrect(:, :), brhs(:, :), work(:, :) ! (nvar,nz) real(fltp), pointer :: x1(:) ! =(nvar,nz) integer, pointer :: ipiv1(:) ! =(nvar,nz) @@ -82,8 +56,7 @@ subroutine test_block(which_decsol_option, for_release) integer, pointer :: ipar_decsol(:) ! (lid) real(fltp) :: time_factor, time_solve, time_refine, time_dealloc - integer :: i, j, k, ierr, lid, lrd, nvar, nz, rep - logical :: use_given_weights + integer :: ierr, lid, lrd, nvar, nz character(len=255) :: fname, which_decsol_str include 'formats' @@ -98,14 +71,11 @@ subroutine test_block(which_decsol_option, for_release) else fname = 'block_tri_12.data' end if - nrep = 1 time_factor = 0; time_solve = 0; time_refine = 0; time_dealloc = 0 call read_testfile(fname) - !call xread_testfile(fname) -#ifdef DBLE if (which_decsol_option == bcyclic_dble) then write (*, *) 'bcyclic_dble' call bcyclic_dble_work_sizes(nvar, nz, lrd, lid) @@ -113,12 +83,10 @@ subroutine test_block(which_decsol_option, for_release) write (*, *) 'bad value for which_decsol_option in test_block' call mesa_error(__FILE__, __LINE__) end if -#endif allocate ( & rpar_decsol(lrd), ipar_decsol(lid), x1(nvar*nz), xcorrect(nvar, nz), & - brhs(nvar, nz), ipiv1(nvar*nz), work(nvar, nz), & - l_init(nvar, nvar, nz), d_init(nvar, nvar, nz), u_init(nvar, nvar, nz), stat=ierr) + brhs(nvar, nz), ipiv1(nvar*nz), work(nvar, nz), stat=ierr) if (ierr /= 0) then write (*, *) 'failed in alloc' call mesa_error(__FILE__, __LINE__) @@ -126,87 +94,34 @@ subroutine test_block(which_decsol_option, for_release) ipiv(1:nvar, 1:nz) => ipiv1(1:nvar*nz) x(1:nvar, 1:nz) => x1(1:nvar*nz) - do k = 1, nz - do i = 1, nvar - do j = 1, nvar - l_init(j, i, k) = lblk(j, i, k) - d_init(j, i, k) = dblk(j, i, k) - u_init(j, i, k) = ublk(j, i, k) - end do - end do - end do - call set_xcorrect call set_brhs(lblk, dblk, ublk) - if (.false.) then ! check data - write (*, *) trim(fname) - write (*, 3) 'nvar nz', nvar, nz - do k = 1, nz - do i = 1, nvar - write (*, 3) 'sol-xcorrect', i, k, brhs(i, k), xcorrect(i, k) - do j = 1, nvar - write (*, 4) 'u-d-l', i, j, k, & - ublk(i, j, k), dblk(i, j, k), lblk(i, j, k) - end do - end do - end do - stop 'test_block_tridiagonal' + if (which_decsol_option == bcyclic_dble) then + call solve_blocks(bcyclic_dble_decsolblk) + else + write (*, *) 'missing case for which_decsol_option', which_decsol_option + call mesa_error(__FILE__, __LINE__) end if - do rep = 1, nrep - - do k = 1, nz - do i = 1, nvar - do j = 1, nvar - lblk(j, i, k) = l_init(j, i, k) - dblk(j, i, k) = d_init(j, i, k) - ublk(j, i, k) = u_init(j, i, k) - end do - end do - end do - - use_given_weights = (rep > 1) - -#ifdef DBLE - if (which_decsol_option == bcyclic_dble) then - call solve_blocks( & - use_given_weights, bcyclic_dble_decsolblk, null_decsolblk_quad) - else - write (*, *) 'missing case for which_decsol_option', which_decsol_option - call mesa_error(__FILE__, __LINE__) - end if - -#endif - - end do - call check_x -#ifdef DBLE if (which_decsol_option == bcyclic_dble) then write (*, *) 'done bcyclic_dble' else if (which_decsol_option == block_thomas_dble) then write (*, *) 'done block_thomas_dble' end if -#else - if (which_decsol_option == block_thomas_quad) then - write (*, *) 'done block_thomas_quad' - end if -#endif write (*, *) deallocate (rpar_decsol, ipar_decsol, x1, xcorrect, work, & - brhs, ipiv1, lblk1, dblk1, ublk1, l_init, d_init, u_init) + brhs, ipiv1, lblk1, dblk1, ublk1) contains - subroutine solve_blocks(use_given_weights, decsolblk, decsolblk_quad) - logical, intent(in) :: use_given_weights + subroutine solve_blocks(decsolblk) interface include 'mtx_decsolblk_dble.dek' - include 'mtx_decsolblk_quad.dek' end interface integer :: iop, rep @@ -215,13 +130,8 @@ subroutine solve_blocks(use_given_weights, decsolblk, decsolblk_quad) include 'formats' iop = 0 ! factor A -#ifdef DBLE call decsolblk( & iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr) -#else - call decsolblk_quad( & - iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr) -#endif if (ierr /= 0) then write (*, *) 'decsolblk failed for factor' call mesa_error(__FILE__, __LINE__) @@ -236,13 +146,8 @@ subroutine solve_blocks(use_given_weights, decsolblk, decsolblk_quad) x(j, k) = brhs(j, k) end do end do -#ifdef DBLE call decsolblk( & iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr) -#else - call decsolblk_quad( & - iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr) -#endif if (ierr /= 0) then write (*, *) 'decsolblk failed for solve' call mesa_error(__FILE__, __LINE__) @@ -251,13 +156,8 @@ subroutine solve_blocks(use_given_weights, decsolblk, decsolblk_quad) end do iop = 2 ! deallocate -#ifdef DBLE call decsolblk( & iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr) -#else - call decsolblk_quad( & - iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr) -#endif if (ierr /= 0) then write (*, *) 'decsolblk failed for deallocate' call mesa_error(__FILE__, __LINE__) @@ -275,11 +175,9 @@ subroutine read_testfile(fname) write (*, *) 'failed to open '//trim(fname) call mesa_error(__FILE__, __LINE__) end if -#ifdef DBLE + call mtx_read_block_tridiagonal(iounit, nvar, nz, lblk1, dblk1, ublk1, ierr) -#else - call mtx_read_quad_block_tridiagonal(iounit, nvar, nz, lblk1, dblk1, ublk1, ierr) -#endif + if (ierr /= 0) then write (*, *) 'failed to read '//trim(fname) call mesa_error(__FILE__, __LINE__) @@ -301,11 +199,9 @@ subroutine set_brhs(lblk, dblk, ublk) integer :: k, j include 'formats' ! set brhs = A*xcorrect -#ifdef DBLE + call block_dble_mv(nvar, nz, lblk, dblk, ublk, xcorrect, brhs) -#else - call block_quad_mv(lblk, dblk, ublk, xcorrect, brhs) -#endif + return do k = 1, 2 !nz do j = 1, nvar @@ -369,9 +265,4 @@ subroutine set_xcorrect end subroutine set_xcorrect end subroutine test_block - -#ifdef DBLE - end module test_block_tri_dble -#else -end module test_block_tri_quad -#endif +end module test_block_tri_dble diff --git a/mtx/test/src/test_mtx.f90 b/mtx/test/src/test_mtx.f90 index 66173b251..42c8b8e6a 100644 --- a/mtx/test/src/test_mtx.f90 +++ b/mtx/test/src/test_mtx.f90 @@ -29,7 +29,6 @@ program test_mtx use test_square use test_block_tri_dble, only: do_test_block_tri_dble - use test_block_tri_quad, only: do_test_block_tri_quad use utils_lib, only: mesa_error @@ -49,7 +48,6 @@ program test_mtx call do_test_square call do_test_block_tri_dble - call do_test_block_tri_quad call test_format_conversion end program From e75260fc5e79a848bc2e8d37b4546ea0044e95d0 Mon Sep 17 00:00:00 2001 From: Vincent Vanlaer Date: Thu, 8 Aug 2024 23:10:06 -0400 Subject: [PATCH 3/3] num, star: cleanup after mtx quad precision removal --- num/public/num_def.f | 1 - star/private/star_utils.f90 | 6 ------ 2 files changed, 7 deletions(-) diff --git a/num/public/num_def.f b/num/public/num_def.f index d7fcc100a..d56e5e64f 100644 --- a/num/public/num_def.f +++ b/num/public/num_def.f @@ -56,7 +56,6 @@ module num_def integer, parameter :: square_matrix_type = 1 integer, parameter :: banded_matrix_type = 2 integer, parameter :: block_tridiag_dble_matrix_type = 3 - integer, parameter :: block_tridiag_quad_matrix_type = 4 ! parameter indices in newton iwork array diff --git a/star/private/star_utils.f90 b/star/private/star_utils.f90 index 346260621..b305c28f2 100644 --- a/star/private/star_utils.f90 +++ b/star/private/star_utils.f90 @@ -2845,8 +2845,6 @@ end subroutine show_matrix ! e00(i,j,k) is partial of equ(i,k) wrt var(j,k) subroutine e00(s,i,j,k,nvar,v) - use num_def, only: & - block_tridiag_dble_matrix_type, block_tridiag_quad_matrix_type type (star_info), pointer :: s integer, intent(in) :: i, j, k, nvar real(dp), intent(in) :: v @@ -2903,8 +2901,6 @@ end subroutine e00 ! em1(i,j,k) is partial of equ(i,k) wrt var(j,k-1) subroutine em1(s,i,j,k,nvar,v) - use num_def, only: & - block_tridiag_dble_matrix_type, block_tridiag_quad_matrix_type type (star_info), pointer :: s integer, intent(in) :: i, j, k, nvar real(dp), intent(in) :: v @@ -2958,8 +2954,6 @@ end subroutine em1 ! ep1(i,j,k) is partial of equ(i,k) wrt var(j,k+1) subroutine ep1(s,i,j,k,nvar,v) - use num_def, only: & - block_tridiag_dble_matrix_type, block_tridiag_quad_matrix_type type (star_info), pointer :: s integer, intent(in) :: i, j, k, nvar real(dp), intent(in) :: v