diff --git a/BLAS/src/dtrsm_b.f b/BLAS/src/dtrsm_b.f new file mode 100644 index 0000000..79bbe86 --- /dev/null +++ b/BLAS/src/dtrsm_b.f @@ -0,0 +1,905 @@ +C Generated by TAPENADE (INRIA, Ecuador team) +C Tapenade 3.16 (develop) - 14 Mar 2025 11:54 +C +C Differentiation of dtrsm in reverse (adjoint) mode: +C gradient of useful results: b +C with respect to varying inputs: alpha a b +C RW status of diff variables: alpha:out a:out b:in-out +C> \brief \b DTRSM +C +C =========== DOCUMENTATION =========== +C +C Online html documentation available at +C http://www.netlib.org/lapack/explore-html/ +C +C Definition: +C =========== +C +C SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +C +C .. Scalar Arguments .. +C DOUBLE PRECISION ALPHA +C INTEGER LDA,LDB,M,N +C CHARACTER DIAG,SIDE,TRANSA,UPLO +C .. +C .. Array Arguments .. +C DOUBLE PRECISION A(LDA,*),B(LDB,*) +C .. +C +C +C> \par Purpose: +C ============= +C> +C> \verbatim +C> +C> DTRSM solves one of the matrix equations +C> +C> op( A )*X = alpha*B, or X*op( A ) = alpha*B, +C> +C> where alpha is a scalar, X and B are m by n matrices, A is a unit, or +C> non-unit, upper or lower triangular matrix and op( A ) is one of +C> +C> op( A ) = A or op( A ) = A**T. +C> +C> The matrix X is overwritten on B. +C> \endverbatim +C +C Arguments: +C ========== +C +C> \param[in] SIDE +C> \verbatim +C> SIDE is CHARACTER*1 +C> On entry, SIDE specifies whether op( A ) appears on the left +C> or right of X as follows: +C> +C> SIDE = 'L' or 'l' op( A )*X = alpha*B. +C> +C> SIDE = 'R' or 'r' X*op( A ) = alpha*B. +C> \endverbatim +C> +C> \param[in] UPLO +C> \verbatim +C> UPLO is CHARACTER*1 +C> On entry, UPLO specifies whether the matrix A is an upper or +C> lower triangular matrix as follows: +C> +C> UPLO = 'U' or 'u' A is an upper triangular matrix. +C> +C> UPLO = 'L' or 'l' A is a lower triangular matrix. +C> \endverbatim +C> +C> \param[in] TRANSA +C> \verbatim +C> TRANSA is CHARACTER*1 +C> On entry, TRANSA specifies the form of op( A ) to be used in +C> the matrix multiplication as follows: +C> +C> TRANSA = 'N' or 'n' op( A ) = A. +C> +C> TRANSA = 'T' or 't' op( A ) = A**T. +C> +C> TRANSA = 'C' or 'c' op( A ) = A**T. +C> \endverbatim +C> +C> \param[in] DIAG +C> \verbatim +C> DIAG is CHARACTER*1 +C> On entry, DIAG specifies whether or not A is unit triangular +C> as follows: +C> +C> DIAG = 'U' or 'u' A is assumed to be unit triangular. +C> +C> DIAG = 'N' or 'n' A is not assumed to be unit +C> triangular. +C> \endverbatim +C> +C> \param[in] M +C> \verbatim +C> M is INTEGER +C> On entry, M specifies the number of rows of B. M must be at +C> least zero. +C> \endverbatim +C> +C> \param[in] N +C> \verbatim +C> N is INTEGER +C> On entry, N specifies the number of columns of B. N must be +C> at least zero. +C> \endverbatim +C> +C> \param[in] ALPHA +C> \verbatim +C> ALPHA is DOUBLE PRECISION. +C> On entry, ALPHA specifies the scalar alpha. When alpha is +C> zero then A is not referenced and B need not be set before +C> entry. +C> \endverbatim +C> +C> \param[in] A +C> \verbatim +C> A is DOUBLE PRECISION array, dimension ( LDA, k ), +C> where k is m when SIDE = 'L' or 'l' +C> and k is n when SIDE = 'R' or 'r'. +C> Before entry with UPLO = 'U' or 'u', the leading k by k +C> upper triangular part of the array A must contain the upper +C> triangular matrix and the strictly lower triangular part of +C> A is not referenced. +C> Before entry with UPLO = 'L' or 'l', the leading k by k +C> lower triangular part of the array A must contain the lower +C> triangular matrix and the strictly upper triangular part of +C> A is not referenced. +C> Note that when DIAG = 'U' or 'u', the diagonal elements of +C> A are not referenced either, but are assumed to be unity. +C> \endverbatim +C> +C> \param[in] LDA +C> \verbatim +C> LDA is INTEGER +C> On entry, LDA specifies the first dimension of A as declared +C> in the calling (sub) program. When SIDE = 'L' or 'l' then +C> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +C> then LDA must be at least max( 1, n ). +C> \endverbatim +C> +C> \param[in,out] B +C> \verbatim +C> B is DOUBLE PRECISION array, dimension ( LDB, N ) +C> Before entry, the leading m by n part of the array B must +C> contain the right-hand side matrix B, and on exit is +C> overwritten by the solution matrix X. +C> \endverbatim +C> +C> \param[in] LDB +C> \verbatim +C> LDB is INTEGER +C> On entry, LDB specifies the first dimension of B as declared +C> in the calling (sub) program. LDB must be at least +C> max( 1, m ). +C> \endverbatim +C +C Authors: +C ======== +C +C> \author Univ. of Tennessee +C> \author Univ. of California Berkeley +C> \author Univ. of Colorado Denver +C> \author NAG Ltd. +C +C> \ingroup trsm +C +C> \par Further Details: +C ===================== +C> +C> \verbatim +C> +C> Level 3 Blas routine. +C> +C> +C> -- Written on 8-February-1989. +C> Jack Dongarra, Argonne National Laboratory. +C> Iain Duff, AERE Harwell. +C> Jeremy Du Croz, Numerical Algorithms Group Ltd. +C> Sven Hammarling, Numerical Algorithms Group Ltd. +C> \endverbatim +C> +C ===================================================================== +C SUBROUTINE DTRSM_B(side, uplo, transa, diag, m, n, alpha, alphab, +C + a, ab, lda, b, bb, ldb) +C IMPLICIT NONE +C INCLUDE 'DIFFSIZES.inc' +CC Hint: ISIZE2OFa should be the size of dimension 2 of array a +CC +CC -- Reference BLAS level3 routine -- +CC -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +CC -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +CC +CC .. Scalar Arguments .. +C DOUBLE PRECISION alpha +C DOUBLE PRECISION alphab +C INTEGER lda, ldb, m, n +C CHARACTER diag, side, transa, uplo +CC .. +CC .. Array Arguments .. +C DOUBLE PRECISION a(lda, *), b(ldb, *) +C DOUBLE PRECISION ab(lda, *), bb(ldb, *) +C EXTERNAL LSAME +CC .. +CC +CC ===================================================================== +CC +CC .. External Functions .. +C INTEGER get_ISIZE2OFA +C EXTERNAL get_ISIZE2OFA +C LOGICAL LSAME +CC .. +CC .. External Subroutines .. +C EXTERNAL XERBLA, check_ISIZE2OFA_initialized +CC .. +CC .. Intrinsic Functions .. +C INTRINSIC MAX +CC .. +CC .. Local Scalars .. +C DOUBLE PRECISION temp +C DOUBLE PRECISION tempb +C INTEGER i, info, j, k, nrowa +C LOGICAL lside, nounit, upper +C INTEGER ISIZE2OFA +CC .. +CC .. Parameters .. +C DOUBLE PRECISION zero +C PARAMETER (zero=0.0d+0) +C INTEGER max1 +C INTEGER max2 +C DOUBLE PRECISION tempb0 +C DOUBLE PRECISION tmp +C DOUBLE PRECISION tmpb +C DOUBLE PRECISION tmp0 +C DOUBLE PRECISION tmpb0 +C DOUBLE PRECISION tmp1 +C DOUBLE PRECISION tmpb1 +C DOUBLE PRECISION tmp2 +C DOUBLE PRECISION tmpb2 +C DOUBLE PRECISION tmp3 +C DOUBLE PRECISION tmpb3 +C DOUBLE PRECISION tmp4 +C DOUBLE PRECISION tmpb4 +C INTEGER ad_to +C INTEGER*4 branch +C INTEGER ad_from +C INTEGER ad_to0 +C INTEGER ad_from0 +C INTEGER ad_to1 +C INTEGER ad_from1 +C INTEGER ad_to2 +C INTEGER ad_from2 +C INTEGER ii2 +C INTEGER ii1 +CC .. +CC +CC Test the input parameters. +CC +C CALL check_ISIZE2OFA_initialized() +C ISIZE2OFA = get_ISIZE2OFA() +C lside = LSAME(side, 'L') +C IF (lside) THEN +C nrowa = m +C ELSE +C nrowa = n +C END IF +C nounit = LSAME(diag, 'N') +C upper = LSAME(uplo, 'U') +CC +C info = 0 +C IF (.NOT.lside .AND. (.NOT.LSAME(side, 'R'))) THEN +C CALL PUSHCONTROL3B(0) +C info = 1 +C ELSE IF (.NOT.upper .AND. (.NOT.LSAME(uplo, 'L'))) THEN +C CALL PUSHCONTROL3B(1) +C info = 2 +C ELSE IF (.NOT.LSAME(transa, 'N') .AND. (.NOT.LSAME(transa, 'T')) +C + .AND. (.NOT.LSAME(transa, 'C'))) THEN +C CALL PUSHCONTROL3B(2) +C info = 3 +C ELSE IF (.NOT.LSAME(diag, 'U') .AND. (.NOT.LSAME(diag, 'N'))) THEN +C CALL PUSHCONTROL3B(3) +C info = 4 +C ELSE IF (m .LT. 0) THEN +C CALL PUSHCONTROL3B(4) +C info = 5 +C ELSE IF (n .LT. 0) THEN +C CALL PUSHCONTROL3B(5) +C info = 6 +C ELSE +C IF (1 .LT. nrowa) THEN +C max1 = nrowa +C ELSE +C max1 = 1 +C END IF +C IF (lda .LT. max1) THEN +C CALL PUSHCONTROL3B(6) +C info = 9 +C ELSE +C IF (1 .LT. m) THEN +C max2 = m +C ELSE +C max2 = 1 +C END IF +C IF (ldb .LT. max2) THEN +C CALL PUSHCONTROL3B(7) +C info = 11 +C ELSE +C CALL PUSHCONTROL3B(7) +C END IF +C END IF +C END IF +C IF (info .EQ. 0) THEN +CC +CC Quick return if possible. +CC +C IF (m .EQ. 0 .OR. n .EQ. 0) THEN +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ELSE IF (alpha .EQ. zero) THEN +CC +CC And when alpha.eq.zero. +CC +C DO j=n,1,-1 +C DO i=m,1,-1 +C bb(i, j) = 0.D0 +C ENDDO +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ELSE IF (lside) THEN +CC +CC Start the operations. +CC +C IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*inv( A )*B. +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=m,1,-1 +C IF (nounit) THEN +C CALL PUSHREAL8(b(k, j)) +C b(k, j) = b(k, j)/a(k, k) +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C DO i=1,k-1 +C tmp = b(i, j) - b(k, j)*a(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp +C ENDDO +C CALL PUSHINTEGER4(i - 1) +C ENDDO +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO k=1,m,1 +C CALL POPINTEGER4(ad_to) +C DO i=ad_to,1,-1 +C CALL POPREAL8(b(i, j)) +C tmpb = bb(i, j) +C bb(i, j) = tmpb +C bb(k, j) = bb(k, j) - a(i, k)*tmpb +C ab(i, k) = ab(i, k) - b(k, j)*tmpb +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C CALL POPREAL8(b(k, j)) +C tempb0 = bb(k, j)/a(k, k) +C bb(k, j) = tempb0 +C ab(k, k) = ab(k, k) - b(k, j)*tempb0/a(k, k) +C END IF +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C alphab = alphab + b(i, j)*bb(i, j) +C bb(i, j) = alpha*bb(i, j) +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,m +C IF (nounit) THEN +C CALL PUSHREAL8(b(k, j)) +C b(k, j) = b(k, j)/a(k, k) +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C ad_from = k + 1 +C DO i=ad_from,m +C tmp0 = b(i, j) - b(k, j)*a(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp0 +C ENDDO +C CALL PUSHINTEGER4(ad_from) +C ENDDO +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO k=m,1,-1 +C CALL POPINTEGER4(ad_from) +C DO i=m,ad_from,-1 +C CALL POPREAL8(b(i, j)) +C tmpb0 = bb(i, j) +C bb(i, j) = tmpb0 +C bb(k, j) = bb(k, j) - a(i, k)*tmpb0 +C ab(i, k) = ab(i, k) - b(k, j)*tmpb0 +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C CALL POPREAL8(b(k, j)) +C tempb0 = bb(k, j)/a(k, k) +C bb(k, j) = tempb0 +C ab(k, k) = ab(k, k) - b(k, j)*tempb0/a(k, k) +C END IF +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C alphab = alphab + b(i, j)*bb(i, j) +C bb(i, j) = alpha*bb(i, j) +C ENDDO +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*inv( A**T )*B. +CC +C DO j=1,n +C DO i=1,m +C temp = alpha*b(i, j) +C DO k=1,i-1 +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C CALL PUSHINTEGER4(k - 1) +C IF (nounit) THEN +C CALL PUSHREAL8(temp) +C temp = temp/a(i, i) +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = temp +C ENDDO +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C tempb = bb(i, j) +C bb(i, j) = 0.D0 +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C CALL POPREAL8(temp) +C tempb0 = tempb/a(i, i) +C tempb = tempb0 +C ab(i, i) = ab(i, i) - temp*tempb0/a(i, i) +C END IF +C CALL POPINTEGER4(ad_to0) +C DO k=ad_to0,1,-1 +C ab(k, i) = ab(k, i) - b(k, j)*tempb +C bb(k, j) = bb(k, j) - a(k, i)*tempb +C ENDDO +C alphab = alphab + b(i, j)*tempb +C bb(i, j) = bb(i, j) + alpha*tempb +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=m,1,-1 +C temp = alpha*b(i, j) +C ad_from0 = i + 1 +C DO k=ad_from0,m +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C CALL PUSHINTEGER4(ad_from0) +C IF (nounit) THEN +C CALL PUSHREAL8(temp) +C temp = temp/a(i, i) +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = temp +C ENDDO +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO i=1,m,1 +C CALL POPREAL8(b(i, j)) +C tempb = bb(i, j) +C bb(i, j) = 0.D0 +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C CALL POPREAL8(temp) +C tempb0 = tempb/a(i, i) +C tempb = tempb0 +C ab(i, i) = ab(i, i) - temp*tempb0/a(i, i) +C END IF +C CALL POPINTEGER4(ad_from0) +C DO k=m,ad_from0,-1 +C ab(k, i) = ab(k, i) - b(k, j)*tempb +C bb(k, j) = bb(k, j) - a(k, i)*tempb +C ENDDO +C alphab = alphab + b(i, j)*tempb +C bb(i, j) = bb(i, j) + alpha*tempb +C ENDDO +C ENDDO +C END IF +C ELSE IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*B*inv( A ). +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,j-1 +C DO i=1,m +C tmp1 = b(i, j) - a(k, j)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp1 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(k - 1) +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = b(i, j)/a(j, j) +C ENDDO +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO j=n,1,-1 +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C tempb0 = bb(i, j)/a(j, j) +C bb(i, j) = tempb0 +C ab(j, j) = ab(j, j) - b(i, j)*tempb0/a(j, j) +C ENDDO +C END IF +C CALL POPINTEGER4(ad_to1) +C DO k=ad_to1,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C tmpb1 = bb(i, j) +C bb(i, j) = tmpb1 +C ab(k, j) = ab(k, j) - b(i, k)*tmpb1 +C bb(i, k) = bb(i, k) - a(k, j)*tmpb1 +C ENDDO +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C alphab = alphab + b(i, j)*bb(i, j) +C bb(i, j) = alpha*bb(i, j) +C ENDDO +C ENDDO +C ELSE +C DO j=n,1,-1 +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C ad_from1 = j + 1 +C DO k=ad_from1,n +C DO i=1,m +C tmp2 = b(i, j) - a(k, j)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp2 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(ad_from1) +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = b(i, j)/a(j, j) +C ENDDO +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO j=1,n,1 +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C tempb0 = bb(i, j)/a(j, j) +C bb(i, j) = tempb0 +C ab(j, j) = ab(j, j) - b(i, j)*tempb0/a(j, j) +C ENDDO +C END IF +C CALL POPINTEGER4(ad_from1) +C DO k=n,ad_from1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C tmpb2 = bb(i, j) +C bb(i, j) = tmpb2 +C ab(k, j) = ab(k, j) - b(i, k)*tmpb2 +C bb(i, k) = bb(i, k) - a(k, j)*tmpb2 +C ENDDO +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C alphab = alphab + b(i, j)*bb(i, j) +C bb(i, j) = alpha*bb(i, j) +C ENDDO +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*B*inv( A**T ). +CC +C DO k=n,1,-1 +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = b(i, k)/a(k, k) +C ENDDO +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C DO j=1,k-1 +C DO i=1,m +C tmp3 = b(i, j) - a(j, k)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp3 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(j - 1) +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO k=1,n,1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C alphab = alphab + b(i, k)*bb(i, k) +C bb(i, k) = alpha*bb(i, k) +C ENDDO +C CALL POPINTEGER4(ad_to2) +C DO j=ad_to2,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C tmpb3 = bb(i, j) +C bb(i, j) = tmpb3 +C ab(j, k) = ab(j, k) - b(i, k)*tmpb3 +C bb(i, k) = bb(i, k) - a(j, k)*tmpb3 +C ENDDO +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C tempb0 = bb(i, k)/a(k, k) +C bb(i, k) = tempb0 +C ab(k, k) = ab(k, k) - b(i, k)*tempb0/a(k, k) +C ENDDO +C END IF +C ENDDO +C ELSE +C DO k=1,n +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = b(i, k)/a(k, k) +C ENDDO +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C ad_from2 = k + 1 +C DO j=ad_from2,n +C DO i=1,m +C tmp4 = b(i, j) - a(j, k)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp4 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(ad_from2) +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C alphab = 0.D0 +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C ab(ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C DO k=n,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C alphab = alphab + b(i, k)*bb(i, k) +C bb(i, k) = alpha*bb(i, k) +C ENDDO +C CALL POPINTEGER4(ad_from2) +C DO j=n,ad_from2,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C tmpb4 = bb(i, j) +C bb(i, j) = tmpb4 +C ab(j, k) = ab(j, k) - b(i, k)*tmpb4 +C bb(i, k) = bb(i, k) - a(j, k)*tmpb4 +C ENDDO +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C tempb0 = bb(i, k)/a(k, k) +C bb(i, k) = tempb0 +C ab(k, k) = ab(k, k) - b(i, k)*tempb0/a(k, k) +C ENDDO +C END IF +C ENDDO +C END IF +C END IF +C CALL POPCONTROL3B(branch) +C END + + SUBROUTINE DTRSM_B(SIDE, UPLO, TRANSA, DIAG, M, N, + + ALPHA, ALPHAB, A, AB, LDA, B, BB, LDB) +C +C Reverse-mode (adjoint) derivative of DTRSM via the Giles/black-box +C method: differentiates op(A)*X = alpha*B directly, using only the +C original, undifferentiated DTRSM -- never the substitution algorithm. +C +C CONVENTION (matches Tapenade's own generated dtrsm_b.f): +C B : untouched, original forward-call input. +C BB : IN = seed adjoint dF/dX; OUT = adjoint dF/dB (in place). +C ALPHAB, AB : OUTPUT ONLY -- zeroed and filled here, not accumulated +C into whatever the caller passed in. +C + IMPLICIT NONE + CHARACTER SIDE, UPLO, TRANSA, DIAG + INTEGER M, N, LDA, LDB + DOUBLE PRECISION ALPHA, ALPHAB + DOUBLE PRECISION A(LDA,*), AB(LDA,*) + DOUBLE PRECISION B(LDB,*), BB(LDB,*) + + DOUBLE PRECISION X(M,N) + CHARACTER TRANSA_T + LOGICAL LEFT, UPPER, UNIT, NOTRANS + INTEGER I, J, K + + LEFT = (SIDE.EQ.'L' .OR. SIDE.EQ.'l') + UPPER = (UPLO.EQ.'U' .OR. UPLO.EQ.'u') + UNIT = (DIAG.EQ.'U' .OR. DIAG.EQ.'u') + NOTRANS = (TRANSA.EQ.'N' .OR. TRANSA.EQ.'n') + IF (NOTRANS) THEN + TRANSA_T = 'T' + ELSE + TRANSA_T = 'N' + END IF + + IF (M.EQ.0 .OR. N.EQ.0) THEN + ALPHAB = 0.0D0 + RETURN + END IF + +C Step 1: Yb solves op(A)^T*Yb = Xb (SIDE='L') or Yb*op(A)^T = Xb (SIDE='R'). +C Overwrites BB in place: held the seed on entry, holds Yb after. + CALL DTRSM(SIDE, UPLO, TRANSA_T, DIAG, M, N, 1.0D0, A, LDA, BB, LDB) + +C Step 2: recompute X from the ORIGINAL, untouched B + DO J = 1, N + DO I = 1, M + X(I,J) = B(I,J) + END DO + END DO + CALL DTRSM(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, X, M) + +C Step 3: alphab = (ORIGINAL B, not X) + ALPHAB = 0.0D0 + DO J = 1, N + DO I = 1, M + ALPHAB = ALPHAB + BB(I,J)*B(I,J) + END DO + END DO + +C Step 4: Ab, restricted to the triangle DTRSM actually reads + DO J = 1, LDA + DO I = 1, LDA + AB(I,J) = 0.0D0 + END DO + END DO + + IF (LEFT) THEN + DO J = 1, M + DO I = 1, M + IF ((UPPER.AND.I.LE.J) .OR. (.NOT.UPPER.AND.I.GE.J)) THEN + IF (.NOT.(UNIT.AND.I.EQ.J)) THEN + IF (NOTRANS) THEN + DO K = 1, N + AB(I,J) = AB(I,J) - BB(I,K)*X(J,K) + END DO + ELSE + DO K = 1, N + AB(I,J) = AB(I,J) - X(I,K)*BB(J,K) + END DO + END IF + END IF + END IF + END DO + END DO + ELSE + DO J = 1, N + DO I = 1, N + IF ((UPPER.AND.I.LE.J) .OR. (.NOT.UPPER.AND.I.GE.J)) THEN + IF (.NOT.(UNIT.AND.I.EQ.J)) THEN + IF (NOTRANS) THEN + DO K = 1, M + AB(I,J) = AB(I,J) - X(K,I)*BB(K,J) + END DO + ELSE + DO K = 1, M + AB(I,J) = AB(I,J) - BB(K,I)*X(K,J) + END DO + END IF + END IF + END IF + END DO + END DO + END IF + +C Step 5: Bb = alpha * Yb (BB currently holds Yb from Step 1) + DO J = 1, N + DO I = 1, M + BB(I,J) = ALPHA*BB(I,J) + END DO + END DO + + RETURN + END SUBROUTINE DTRSM_B diff --git a/BLAS/src/dtrsm_bv.f b/BLAS/src/dtrsm_bv.f new file mode 100644 index 0000000..694e6b2 --- /dev/null +++ b/BLAS/src/dtrsm_bv.f @@ -0,0 +1,1054 @@ +C Generated by TAPENADE (INRIA, Ecuador team) +C Tapenade 3.16 (develop) - 14 Mar 2025 11:54 +C +C Differentiation of dtrsm in reverse (adjoint) mode (with options multiDirectional): +C gradient of useful results: b +C with respect to varying inputs: alpha a b +C RW status of diff variables: alpha:out a:out b:in-out +C> \brief \b DTRSM +C +C =========== DOCUMENTATION =========== +C +C Online html documentation available at +C http://www.netlib.org/lapack/explore-html/ +C +C Definition: +C =========== +C +C SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +C +C .. Scalar Arguments .. +C DOUBLE PRECISION ALPHA +C INTEGER LDA,LDB,M,N +C CHARACTER DIAG,SIDE,TRANSA,UPLO +C .. +C .. Array Arguments .. +C DOUBLE PRECISION A(LDA,*),B(LDB,*) +C .. +C +C +C> \par Purpose: +C ============= +C> +C> \verbatim +C> +C> DTRSM solves one of the matrix equations +C> +C> op( A )*X = alpha*B, or X*op( A ) = alpha*B, +C> +C> where alpha is a scalar, X and B are m by n matrices, A is a unit, or +C> non-unit, upper or lower triangular matrix and op( A ) is one of +C> +C> op( A ) = A or op( A ) = A**T. +C> +C> The matrix X is overwritten on B. +C> \endverbatim +C +C Arguments: +C ========== +C +C> \param[in] SIDE +C> \verbatim +C> SIDE is CHARACTER*1 +C> On entry, SIDE specifies whether op( A ) appears on the left +C> or right of X as follows: +C> +C> SIDE = 'L' or 'l' op( A )*X = alpha*B. +C> +C> SIDE = 'R' or 'r' X*op( A ) = alpha*B. +C> \endverbatim +C> +C> \param[in] UPLO +C> \verbatim +C> UPLO is CHARACTER*1 +C> On entry, UPLO specifies whether the matrix A is an upper or +C> lower triangular matrix as follows: +C> +C> UPLO = 'U' or 'u' A is an upper triangular matrix. +C> +C> UPLO = 'L' or 'l' A is a lower triangular matrix. +C> \endverbatim +C> +C> \param[in] TRANSA +C> \verbatim +C> TRANSA is CHARACTER*1 +C> On entry, TRANSA specifies the form of op( A ) to be used in +C> the matrix multiplication as follows: +C> +C> TRANSA = 'N' or 'n' op( A ) = A. +C> +C> TRANSA = 'T' or 't' op( A ) = A**T. +C> +C> TRANSA = 'C' or 'c' op( A ) = A**T. +C> \endverbatim +C> +C> \param[in] DIAG +C> \verbatim +C> DIAG is CHARACTER*1 +C> On entry, DIAG specifies whether or not A is unit triangular +C> as follows: +C> +C> DIAG = 'U' or 'u' A is assumed to be unit triangular. +C> +C> DIAG = 'N' or 'n' A is not assumed to be unit +C> triangular. +C> \endverbatim +C> +C> \param[in] M +C> \verbatim +C> M is INTEGER +C> On entry, M specifies the number of rows of B. M must be at +C> least zero. +C> \endverbatim +C> +C> \param[in] N +C> \verbatim +C> N is INTEGER +C> On entry, N specifies the number of columns of B. N must be +C> at least zero. +C> \endverbatim +C> +C> \param[in] ALPHA +C> \verbatim +C> ALPHA is DOUBLE PRECISION. +C> On entry, ALPHA specifies the scalar alpha. When alpha is +C> zero then A is not referenced and B need not be set before +C> entry. +C> \endverbatim +C> +C> \param[in] A +C> \verbatim +C> A is DOUBLE PRECISION array, dimension ( LDA, k ), +C> where k is m when SIDE = 'L' or 'l' +C> and k is n when SIDE = 'R' or 'r'. +C> Before entry with UPLO = 'U' or 'u', the leading k by k +C> upper triangular part of the array A must contain the upper +C> triangular matrix and the strictly lower triangular part of +C> A is not referenced. +C> Before entry with UPLO = 'L' or 'l', the leading k by k +C> lower triangular part of the array A must contain the lower +C> triangular matrix and the strictly upper triangular part of +C> A is not referenced. +C> Note that when DIAG = 'U' or 'u', the diagonal elements of +C> A are not referenced either, but are assumed to be unity. +C> \endverbatim +C> +C> \param[in] LDA +C> \verbatim +C> LDA is INTEGER +C> On entry, LDA specifies the first dimension of A as declared +C> in the calling (sub) program. When SIDE = 'L' or 'l' then +C> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +C> then LDA must be at least max( 1, n ). +C> \endverbatim +C> +C> \param[in,out] B +C> \verbatim +C> B is DOUBLE PRECISION array, dimension ( LDB, N ) +C> Before entry, the leading m by n part of the array B must +C> contain the right-hand side matrix B, and on exit is +C> overwritten by the solution matrix X. +C> \endverbatim +C> +C> \param[in] LDB +C> \verbatim +C> LDB is INTEGER +C> On entry, LDB specifies the first dimension of B as declared +C> in the calling (sub) program. LDB must be at least +C> max( 1, m ). +C> \endverbatim +C +C Authors: +C ======== +C +C> \author Univ. of Tennessee +C> \author Univ. of California Berkeley +C> \author Univ. of Colorado Denver +C> \author NAG Ltd. +C +C> \ingroup trsm +C +C> \par Further Details: +C ===================== +C> +C> \verbatim +C> +C> Level 3 Blas routine. +C> +C> +C> -- Written on 8-February-1989. +C> Jack Dongarra, Argonne National Laboratory. +C> Iain Duff, AERE Harwell. +C> Jeremy Du Croz, Numerical Algorithms Group Ltd. +C> Sven Hammarling, Numerical Algorithms Group Ltd. +C> \endverbatim +C> +C ===================================================================== +C SUBROUTINE DTRSM_BV(side, uplo, transa, diag, m, n, alpha, alphab +C + , a, ab, lda, b, bb, ldb, nbdirs) +C IMPLICIT NONE +C INCLUDE 'DIFFSIZES.inc' +CC Hint: ISIZE2OFa should be the size of dimension 2 of array a +CC Hint: nbdirsmax should be the maximum number of differentiation directions +CC +CC -- Reference BLAS level3 routine -- +CC -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +CC -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +CC +CC .. Scalar Arguments .. +C DOUBLE PRECISION alpha +C DOUBLE PRECISION alphab(nbdirsmax) +C INTEGER lda, ldb, m, n +C CHARACTER diag, side, transa, uplo +CC .. +CC .. Array Arguments .. +C DOUBLE PRECISION a(lda, *), b(ldb, *) +C DOUBLE PRECISION ab(nbdirsmax, lda, *), bb(nbdirsmax, ldb, *) +C EXTERNAL LSAME +CC .. +CC +CC ===================================================================== +CC +CC .. External Functions .. +C INTEGER get_ISIZE2OFA +C EXTERNAL get_ISIZE2OFA +C LOGICAL LSAME +CC .. +CC .. External Subroutines .. +C EXTERNAL XERBLA, check_ISIZE2OFA_initialized +CC .. +CC .. Intrinsic Functions .. +C INTRINSIC MAX +CC .. +CC .. Local Scalars .. +C DOUBLE PRECISION temp +C DOUBLE PRECISION tempb(nbdirsmax) +C INTEGER i, info, j, k, nrowa +C LOGICAL lside, nounit, upper +C INTEGER ISIZE2OFA +CC .. +CC .. Parameters .. +C DOUBLE PRECISION zero +C PARAMETER (zero=0.0d+0) +C INTEGER max1 +C INTEGER max2 +C INTEGER nd +C DOUBLE PRECISION tempb0(nbdirsmax) +C DOUBLE PRECISION tmp +C DOUBLE PRECISION tmpb(nbdirsmax) +C DOUBLE PRECISION tmp0 +C DOUBLE PRECISION tmpb0(nbdirsmax) +C DOUBLE PRECISION tmp1 +C DOUBLE PRECISION tmpb1(nbdirsmax) +C DOUBLE PRECISION tmp2 +C DOUBLE PRECISION tmpb2(nbdirsmax) +C DOUBLE PRECISION tmp3 +C DOUBLE PRECISION tmpb3(nbdirsmax) +C DOUBLE PRECISION tmp4 +C DOUBLE PRECISION tmpb4(nbdirsmax) +C INTEGER ad_to +C INTEGER*4 branch +C INTEGER ad_from +C INTEGER ad_to0 +C INTEGER ad_from0 +C INTEGER ad_to1 +C INTEGER ad_from1 +C INTEGER ad_to2 +C INTEGER ad_from2 +C INTEGER ii2 +C INTEGER ii1 +C INTEGER nbdirs +CC .. +CC +CC Test the input parameters. +CC +CC Check 0 < nbdirs <= nbdirsmax (required by DIFFSIZES.inc) +C CALL check_ISIZE2OFA_initialized() +C ISIZE2OFA = get_ISIZE2OFA() +C IF (nbdirs.LE.0 .OR. nbdirs.GT.nbdirsmax) THEN +C WRITE(*,'(A,I0,A,I0,A)') 'Error: nbdirs=', nbdirs, +C + ' must be in 1..nbdirsmax=', nbdirsmax, '. Stopping.' +C STOP 1 +C END IF +CC +C lside = LSAME(side, 'L') +C IF (lside) THEN +C nrowa = m +C ELSE +C nrowa = n +C END IF +C nounit = LSAME(diag, 'N') +C upper = LSAME(uplo, 'U') +CC +C info = 0 +C IF (.NOT.lside .AND. (.NOT.LSAME(side, 'R'))) THEN +C CALL PUSHCONTROL3B(0) +C info = 1 +C ELSE IF (.NOT.upper .AND. (.NOT.LSAME(uplo, 'L'))) THEN +C CALL PUSHCONTROL3B(1) +C info = 2 +C ELSE IF (.NOT.LSAME(transa, 'N') .AND. (.NOT.LSAME(transa, 'T')) +C + .AND. (.NOT.LSAME(transa, 'C'))) THEN +C CALL PUSHCONTROL3B(2) +C info = 3 +C ELSE IF (.NOT.LSAME(diag, 'U') .AND. (.NOT.LSAME(diag, 'N'))) THEN +C CALL PUSHCONTROL3B(3) +C info = 4 +C ELSE IF (m .LT. 0) THEN +C CALL PUSHCONTROL3B(4) +C info = 5 +C ELSE IF (n .LT. 0) THEN +C CALL PUSHCONTROL3B(5) +C info = 6 +C ELSE +C IF (1 .LT. nrowa) THEN +C max1 = nrowa +C ELSE +C max1 = 1 +C END IF +C IF (lda .LT. max1) THEN +C CALL PUSHCONTROL3B(6) +C info = 9 +C ELSE +C IF (1 .LT. m) THEN +C max2 = m +C ELSE +C max2 = 1 +C END IF +C IF (ldb .LT. max2) THEN +C CALL PUSHCONTROL3B(7) +C info = 11 +C ELSE +C CALL PUSHCONTROL3B(7) +C END IF +C END IF +C END IF +C IF (info .EQ. 0) THEN +CC +CC Quick return if possible. +CC +C IF (m .EQ. 0 .OR. n .EQ. 0) THEN +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C ELSE IF (alpha .EQ. zero) THEN +CC +CC And when alpha.eq.zero. +CC +C DO j=n,1,-1 +C DO i=m,1,-1 +C DO nd=1,nbdirs +C bb(nd, i, j) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C ELSE IF (lside) THEN +CC +CC Start the operations. +CC +C IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*inv( A )*B. +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=m,1,-1 +C IF (nounit) THEN +C CALL PUSHREAL8(b(k, j)) +C b(k, j) = b(k, j)/a(k, k) +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C DO i=1,k-1 +C tmp = b(i, j) - b(k, j)*a(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp +C ENDDO +C CALL PUSHINTEGER4(i - 1) +C ENDDO +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO k=1,m,1 +C CALL POPINTEGER4(ad_to) +C DO i=ad_to,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tmpb(nd) = bb(nd, i, j) +C bb(nd, i, j) = tmpb(nd) +C bb(nd, k, j) = bb(nd, k, j) - a(i, k)*tmpb(nd) +C ab(nd, i, k) = ab(nd, i, k) - b(k, j)*tmpb(nd) +C ENDDO +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C CALL POPREAL8(b(k, j)) +C DO nd=1,nbdirs +C tempb0(nd) = bb(nd, k, j)/a(k, k) +C bb(nd, k, j) = tempb0(nd) +C ab(nd, k, k) = ab(nd, k, k) - b(k, j)*tempb0(nd)/a +C + (k, k) +C ENDDO +C END IF +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, j)*bb(nd, i, j) +C bb(nd, i, j) = alpha*bb(nd, i, j) +C ENDDO +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,m +C IF (nounit) THEN +C CALL PUSHREAL8(b(k, j)) +C b(k, j) = b(k, j)/a(k, k) +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C ad_from = k + 1 +C DO i=ad_from,m +C tmp0 = b(i, j) - b(k, j)*a(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp0 +C ENDDO +C CALL PUSHINTEGER4(ad_from) +C ENDDO +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO k=m,1,-1 +C CALL POPINTEGER4(ad_from) +C DO i=m,ad_from,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tmpb0(nd) = bb(nd, i, j) +C bb(nd, i, j) = tmpb0(nd) +C bb(nd, k, j) = bb(nd, k, j) - a(i, k)*tmpb0(nd) +C ab(nd, i, k) = ab(nd, i, k) - b(k, j)*tmpb0(nd) +C ENDDO +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C CALL POPREAL8(b(k, j)) +C DO nd=1,nbdirs +C tempb0(nd) = bb(nd, k, j)/a(k, k) +C bb(nd, k, j) = tempb0(nd) +C ab(nd, k, k) = ab(nd, k, k) - b(k, j)*tempb0(nd)/a +C + (k, k) +C ENDDO +C END IF +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, j)*bb(nd, i, j) +C bb(nd, i, j) = alpha*bb(nd, i, j) +C ENDDO +C ENDDO +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*inv( A**T )*B. +CC +C DO j=1,n +C DO i=1,m +C temp = alpha*b(i, j) +C DO k=1,i-1 +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C CALL PUSHINTEGER4(k - 1) +C IF (nounit) THEN +C CALL PUSHREAL8(temp) +C temp = temp/a(i, i) +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = temp +C ENDDO +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tempb(nd) = bb(nd, i, j) +C bb(nd, i, j) = 0.D0 +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C CALL POPREAL8(temp) +C DO nd=1,nbdirs +C tempb(nd) = tempb0(nd) +C ab(nd, i, i) = ab(nd, i, i) - temp*tempb0(nd)/a(i, i +C + ) +C ENDDO +C END IF +C CALL POPINTEGER4(ad_to0) +C DO k=ad_to0,1,-1 +C DO nd=1,nbdirs +C ab(nd, k, i) = ab(nd, k, i) - b(k, j)*tempb(nd) +C bb(nd, k, j) = bb(nd, k, j) - a(k, i)*tempb(nd) +C ENDDO +C ENDDO +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, j)*tempb(nd) +C bb(nd, i, j) = bb(nd, i, j) + alpha*tempb(nd) +C ENDDO +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=m,1,-1 +C temp = alpha*b(i, j) +C ad_from0 = i + 1 +C DO k=ad_from0,m +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C CALL PUSHINTEGER4(ad_from0) +C IF (nounit) THEN +C CALL PUSHREAL8(temp) +C temp = temp/a(i, i) +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = temp +C ENDDO +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO j=n,1,-1 +C DO i=1,m,1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tempb(nd) = bb(nd, i, j) +C bb(nd, i, j) = 0.D0 +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C CALL POPREAL8(temp) +C DO nd=1,nbdirs +C tempb(nd) = tempb0(nd) +C ab(nd, i, i) = ab(nd, i, i) - temp*tempb0(nd)/a(i, i +C + ) +C ENDDO +C END IF +C CALL POPINTEGER4(ad_from0) +C DO k=m,ad_from0,-1 +C DO nd=1,nbdirs +C ab(nd, k, i) = ab(nd, k, i) - b(k, j)*tempb(nd) +C bb(nd, k, j) = bb(nd, k, j) - a(k, i)*tempb(nd) +C ENDDO +C ENDDO +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, j)*tempb(nd) +C bb(nd, i, j) = bb(nd, i, j) + alpha*tempb(nd) +C ENDDO +C ENDDO +C ENDDO +C END IF +C ELSE IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*B*inv( A ). +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,j-1 +C DO i=1,m +C tmp1 = b(i, j) - a(k, j)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp1 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(k - 1) +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = b(i, j)/a(j, j) +C ENDDO +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO j=n,1,-1 +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tempb0(nd) = bb(nd, i, j)/a(j, j) +C bb(nd, i, j) = tempb0(nd) +C ab(nd, j, j) = ab(nd, j, j) - b(i, j)*tempb0(nd)/a(j +C + , j) +C ENDDO +C ENDDO +C END IF +C CALL POPINTEGER4(ad_to1) +C DO k=ad_to1,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tmpb1(nd) = bb(nd, i, j) +C bb(nd, i, j) = tmpb1(nd) +C ab(nd, k, j) = ab(nd, k, j) - b(i, k)*tmpb1(nd) +C bb(nd, i, k) = bb(nd, i, k) - a(k, j)*tmpb1(nd) +C ENDDO +C ENDDO +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, j)*bb(nd, i, j) +C bb(nd, i, j) = alpha*bb(nd, i, j) +C ENDDO +C ENDDO +C ENDDO +C ELSE +C DO j=n,1,-1 +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = alpha*b(i, j) +C ENDDO +C ad_from1 = j + 1 +C DO k=ad_from1,n +C DO i=1,m +C tmp2 = b(i, j) - a(k, j)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp2 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(ad_from1) +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = b(i, j)/a(j, j) +C ENDDO +C CALL PUSHCONTROL1B(0) +C ELSE +C CALL PUSHCONTROL1B(1) +C END IF +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO j=1,n,1 +C CALL POPCONTROL1B(branch) +C IF (branch .EQ. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tempb0(nd) = bb(nd, i, j)/a(j, j) +C bb(nd, i, j) = tempb0(nd) +C ab(nd, j, j) = ab(nd, j, j) - b(i, j)*tempb0(nd)/a(j +C + , j) +C ENDDO +C ENDDO +C END IF +C CALL POPINTEGER4(ad_from1) +C DO k=n,ad_from1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tmpb2(nd) = bb(nd, i, j) +C bb(nd, i, j) = tmpb2(nd) +C ab(nd, k, j) = ab(nd, k, j) - b(i, k)*tmpb2(nd) +C bb(nd, i, k) = bb(nd, i, k) - a(k, j)*tmpb2(nd) +C ENDDO +C ENDDO +C ENDDO +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, j)*bb(nd, i, j) +C bb(nd, i, j) = alpha*bb(nd, i, j) +C ENDDO +C ENDDO +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*B*inv( A**T ). +CC +C DO k=n,1,-1 +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = b(i, k)/a(k, k) +C ENDDO +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C DO j=1,k-1 +C DO i=1,m +C tmp3 = b(i, j) - a(j, k)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp3 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(j - 1) +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO k=1,n,1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, k)*bb(nd, i, k) +C bb(nd, i, k) = alpha*bb(nd, i, k) +C ENDDO +C ENDDO +C CALL POPINTEGER4(ad_to2) +C DO j=ad_to2,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tmpb3(nd) = bb(nd, i, j) +C bb(nd, i, j) = tmpb3(nd) +C ab(nd, j, k) = ab(nd, j, k) - b(i, k)*tmpb3(nd) +C bb(nd, i, k) = bb(nd, i, k) - a(j, k)*tmpb3(nd) +C ENDDO +C ENDDO +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C DO nd=1,nbdirs +C tempb0(nd) = bb(nd, i, k)/a(k, k) +C bb(nd, i, k) = tempb0(nd) +C ab(nd, k, k) = ab(nd, k, k) - b(i, k)*tempb0(nd)/a(k, +C + k) +C ENDDO +C ENDDO +C END IF +C ENDDO +C ELSE +C DO k=1,n +C IF (nounit) THEN +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = b(i, k)/a(k, k) +C ENDDO +C CALL PUSHCONTROL1B(1) +C ELSE +C CALL PUSHCONTROL1B(0) +C END IF +C ad_from2 = k + 1 +C DO j=ad_from2,n +C DO i=1,m +C tmp4 = b(i, j) - a(j, k)*b(i, k) +C CALL PUSHREAL8(b(i, j)) +C b(i, j) = tmp4 +C ENDDO +C ENDDO +C CALL PUSHINTEGER4(ad_from2) +C DO i=1,m +C CALL PUSHREAL8(b(i, k)) +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C DO nd=1,nbdirsmax +C alphab(nd) = 0.D0 +C ENDDO +C DO ii1=1,ISIZE2OFa +C DO ii2=1,lda +C DO nd=1,nbdirsmax +C ab(nd, ii2, ii1) = 0.D0 +C ENDDO +C ENDDO +C ENDDO +C DO k=n,1,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C DO nd=1,nbdirs +C alphab(nd) = alphab(nd) + b(i, k)*bb(nd, i, k) +C bb(nd, i, k) = alpha*bb(nd, i, k) +C ENDDO +C ENDDO +C CALL POPINTEGER4(ad_from2) +C DO j=n,ad_from2,-1 +C DO i=m,1,-1 +C CALL POPREAL8(b(i, j)) +C DO nd=1,nbdirs +C tmpb4(nd) = bb(nd, i, j) +C bb(nd, i, j) = tmpb4(nd) +C ab(nd, j, k) = ab(nd, j, k) - b(i, k)*tmpb4(nd) +C bb(nd, i, k) = bb(nd, i, k) - a(j, k)*tmpb4(nd) +C ENDDO +C ENDDO +C ENDDO +C CALL POPCONTROL1B(branch) +C IF (branch .NE. 0) THEN +C DO i=m,1,-1 +C CALL POPREAL8(b(i, k)) +C DO nd=1,nbdirs +C tempb0(nd) = bb(nd, i, k)/a(k, k) +C bb(nd, i, k) = tempb0(nd) +C ab(nd, k, k) = ab(nd, k, k) - b(i, k)*tempb0(nd)/a(k, +C + k) +C ENDDO +C ENDDO +C END IF +C ENDDO +C END IF +C END IF +C CALL POPCONTROL3B(branch) +C END + + SUBROUTINE DTRSM_BV(SIDE, UPLO, TRANSA, DIAG, M, N, + + ALPHA, ALPHAB, A, AB, LDA, B, BB, LDB, NBDIRS) +C +C Vector reverse-mode (adjoint) derivative of DTRSM, black-box/Giles-style. +C Self-contained -- does not call DTRSM_B. AB/BB store direction as the +C FASTEST-varying index (nbdirsmax, LDA/LDB, *): gather each direction into +C a plain work array before use, scatter the result back after. +C +C CONVENTION (matches Tapenade's dtrsm_bv.f): +C B : untouched, original forward-call input. +C BB : IN = seed adjoints (one per direction); OUT = dF/dB. +C ALPHAB, AB : OUTPUT ONLY, zeroed and filled here per direction. +C + IMPLICIT NONE + INCLUDE 'DIFFSIZES.inc' + CHARACTER SIDE, UPLO, TRANSA, DIAG + INTEGER M, N, LDA, LDB, NBDIRS + DOUBLE PRECISION ALPHA + DOUBLE PRECISION ALPHAB(NBDIRSMAX) + DOUBLE PRECISION A(LDA,*), AB(NBDIRSMAX,LDA,*) + DOUBLE PRECISION B(LDB,*), BB(NBDIRSMAX,LDB,*) + + DOUBLE PRECISION X(M,N) + DOUBLE PRECISION BB_SLICE(LDB,N), AB_SLICE(LDA,LDA) + CHARACTER TRANSA_T + LOGICAL LEFT, UPPER, UNIT, NOTRANS + INTEGER ND, I, J, K + + IF (NBDIRS.LE.0 .OR. NBDIRS.GT.NBDIRSMAX) THEN + WRITE(*,'(A,I0,A,I0,A)') 'Error: nbdirs=', NBDIRS, + + ' must be in 1..nbdirsmax=', NBDIRSMAX, '. Stopping.' + STOP 1 + END IF + + LEFT = (SIDE.EQ.'L' .OR. SIDE.EQ.'l') + UPPER = (UPLO.EQ.'U' .OR. UPLO.EQ.'u') + UNIT = (DIAG.EQ.'U' .OR. DIAG.EQ.'u') + NOTRANS = (TRANSA.EQ.'N' .OR. TRANSA.EQ.'n') + IF (NOTRANS) THEN + TRANSA_T = 'T' + ELSE + TRANSA_T = 'N' + END IF + + IF (M.EQ.0 .OR. N.EQ.0) THEN + DO ND = 1, NBDIRS + ALPHAB(ND) = 0.0D0 + END DO + RETURN + END IF + +C X = solution of the original forward call (same for every direction) + DO J = 1, N + DO I = 1, M + X(I,J) = B(I,J) + END DO + END DO + CALL DTRSM(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, X, M) + + DO ND = 1, NBDIRS + +C Gather this direction's seed out of the interleaved BB layout + DO J = 1, N + DO I = 1, M + BB_SLICE(I,J) = BB(ND,I,J) + END DO + END DO + +C Step 1: Yb solves op(A)^T*Yb = Xb (overwrites BB_SLICE in place) + CALL DTRSM(SIDE, UPLO, TRANSA_T, DIAG, M, N, 1.0D0, + + A, LDA, BB_SLICE, LDB) + +C Step 2: alphab(nd) = (ORIGINAL B, not X) + ALPHAB(ND) = 0.0D0 + DO J = 1, N + DO I = 1, M + ALPHAB(ND) = ALPHAB(ND) + BB_SLICE(I,J)*B(I,J) + END DO + END DO + +C Step 3: Ab for this direction, restricted to the triangle read + DO J = 1, LDA + DO I = 1, LDA + AB_SLICE(I,J) = 0.0D0 + END DO + END DO + + IF (LEFT) THEN + DO J = 1, M + DO I = 1, M + IF ((UPPER.AND.I.LE.J).OR.(.NOT.UPPER.AND.I.GE.J)) THEN + IF (.NOT.(UNIT.AND.I.EQ.J)) THEN + IF (NOTRANS) THEN + DO K = 1, N + AB_SLICE(I,J) = AB_SLICE(I,J) + + - BB_SLICE(I,K)*X(J,K) + END DO + ELSE + DO K = 1, N + AB_SLICE(I,J) = AB_SLICE(I,J) + + - X(I,K)*BB_SLICE(J,K) + END DO + END IF + END IF + END IF + END DO + END DO + ELSE + DO J = 1, N + DO I = 1, N + IF ((UPPER.AND.I.LE.J).OR.(.NOT.UPPER.AND.I.GE.J)) THEN + IF (.NOT.(UNIT.AND.I.EQ.J)) THEN + IF (NOTRANS) THEN + DO K = 1, M + AB_SLICE(I,J) = AB_SLICE(I,J) + + - X(K,I)*BB_SLICE(K,J) + END DO + ELSE + DO K = 1, M + AB_SLICE(I,J) = AB_SLICE(I,J) + + - BB_SLICE(K,I)*X(K,J) + END DO + END IF + END IF + END IF + END DO + END DO + END IF + +C Step 4: Bb = alpha * Yb (BB_SLICE currently holds Yb) + DO J = 1, N + DO I = 1, M + BB_SLICE(I,J) = ALPHA*BB_SLICE(I,J) + END DO + END DO + +C Scatter this direction's results back into the interleaved layout + DO J = 1, LDA + DO I = 1, LDA + AB(ND,I,J) = AB_SLICE(I,J) + END DO + END DO + DO J = 1, N + DO I = 1, M + BB(ND,I,J) = BB_SLICE(I,J) + END DO + END DO + + END DO + + RETURN + END SUBROUTINE DTRSM_BV + diff --git a/BLAS/src/dtrsm_d.f b/BLAS/src/dtrsm_d.f new file mode 100644 index 0000000..da09d20 --- /dev/null +++ b/BLAS/src/dtrsm_d.f @@ -0,0 +1,530 @@ +C Generated by TAPENADE (INRIA, Ecuador team) +C Tapenade 3.16 (develop) - 14 Mar 2025 11:54 +C +C Differentiation of dtrsm in forward (tangent) mode: +C variations of useful results: b +C with respect to varying inputs: alpha a b +C RW status of diff variables: alpha:in a:in b:in-out +C> \brief \b DTRSM +C +C =========== DOCUMENTATION =========== +C +C Online html documentation available at +C http://www.netlib.org/lapack/explore-html/ +C +C Definition: +C =========== +C +C SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +C +C .. Scalar Arguments .. +C DOUBLE PRECISION ALPHA +C INTEGER LDA,LDB,M,N +C CHARACTER DIAG,SIDE,TRANSA,UPLO +C .. +C .. Array Arguments .. +C DOUBLE PRECISION A(LDA,*),B(LDB,*) +C .. +C +C +C> \par Purpose: +C ============= +C> +C> \verbatim +C> +C> DTRSM solves one of the matrix equations +C> +C> op( A )*X = alpha*B, or X*op( A ) = alpha*B, +C> +C> where alpha is a scalar, X and B are m by n matrices, A is a unit, or +C> non-unit, upper or lower triangular matrix and op( A ) is one of +C> +C> op( A ) = A or op( A ) = A**T. +C> +C> The matrix X is overwritten on B. +C> \endverbatim +C +C Arguments: +C ========== +C +C> \param[in] SIDE +C> \verbatim +C> SIDE is CHARACTER*1 +C> On entry, SIDE specifies whether op( A ) appears on the left +C> or right of X as follows: +C> +C> SIDE = 'L' or 'l' op( A )*X = alpha*B. +C> +C> SIDE = 'R' or 'r' X*op( A ) = alpha*B. +C> \endverbatim +C> +C> \param[in] UPLO +C> \verbatim +C> UPLO is CHARACTER*1 +C> On entry, UPLO specifies whether the matrix A is an upper or +C> lower triangular matrix as follows: +C> +C> UPLO = 'U' or 'u' A is an upper triangular matrix. +C> +C> UPLO = 'L' or 'l' A is a lower triangular matrix. +C> \endverbatim +C> +C> \param[in] TRANSA +C> \verbatim +C> TRANSA is CHARACTER*1 +C> On entry, TRANSA specifies the form of op( A ) to be used in +C> the matrix multiplication as follows: +C> +C> TRANSA = 'N' or 'n' op( A ) = A. +C> +C> TRANSA = 'T' or 't' op( A ) = A**T. +C> +C> TRANSA = 'C' or 'c' op( A ) = A**T. +C> \endverbatim +C> +C> \param[in] DIAG +C> \verbatim +C> DIAG is CHARACTER*1 +C> On entry, DIAG specifies whether or not A is unit triangular +C> as follows: +C> +C> DIAG = 'U' or 'u' A is assumed to be unit triangular. +C> +C> DIAG = 'N' or 'n' A is not assumed to be unit +C> triangular. +C> \endverbatim +C> +C> \param[in] M +C> \verbatim +C> M is INTEGER +C> On entry, M specifies the number of rows of B. M must be at +C> least zero. +C> \endverbatim +C> +C> \param[in] N +C> \verbatim +C> N is INTEGER +C> On entry, N specifies the number of columns of B. N must be +C> at least zero. +C> \endverbatim +C> +C> \param[in] ALPHA +C> \verbatim +C> ALPHA is DOUBLE PRECISION. +C> On entry, ALPHA specifies the scalar alpha. When alpha is +C> zero then A is not referenced and B need not be set before +C> entry. +C> \endverbatim +C> +C> \param[in] A +C> \verbatim +C> A is DOUBLE PRECISION array, dimension ( LDA, k ), +C> where k is m when SIDE = 'L' or 'l' +C> and k is n when SIDE = 'R' or 'r'. +C> Before entry with UPLO = 'U' or 'u', the leading k by k +C> upper triangular part of the array A must contain the upper +C> triangular matrix and the strictly lower triangular part of +C> A is not referenced. +C> Before entry with UPLO = 'L' or 'l', the leading k by k +C> lower triangular part of the array A must contain the lower +C> triangular matrix and the strictly upper triangular part of +C> A is not referenced. +C> Note that when DIAG = 'U' or 'u', the diagonal elements of +C> A are not referenced either, but are assumed to be unity. +C> \endverbatim +C> +C> \param[in] LDA +C> \verbatim +C> LDA is INTEGER +C> On entry, LDA specifies the first dimension of A as declared +C> in the calling (sub) program. When SIDE = 'L' or 'l' then +C> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +C> then LDA must be at least max( 1, n ). +C> \endverbatim +C> +C> \param[in,out] B +C> \verbatim +C> B is DOUBLE PRECISION array, dimension ( LDB, N ) +C> Before entry, the leading m by n part of the array B must +C> contain the right-hand side matrix B, and on exit is +C> overwritten by the solution matrix X. +C> \endverbatim +C> +C> \param[in] LDB +C> \verbatim +C> LDB is INTEGER +C> On entry, LDB specifies the first dimension of B as declared +C> in the calling (sub) program. LDB must be at least +C> max( 1, m ). +C> \endverbatim +C +C Authors: +C ======== +C +C> \author Univ. of Tennessee +C> \author Univ. of California Berkeley +C> \author Univ. of Colorado Denver +C> \author NAG Ltd. +C +C> \ingroup trsm +C +C> \par Further Details: +C ===================== +C> +C> \verbatim +C> +C> Level 3 Blas routine. +C> +C> +C> -- Written on 8-February-1989. +C> Jack Dongarra, Argonne National Laboratory. +C> Iain Duff, AERE Harwell. +C> Jeremy Du Croz, Numerical Algorithms Group Ltd. +C> Sven Hammarling, Numerical Algorithms Group Ltd. +C> \endverbatim +C> +C ===================================================================== +C SUBROUTINE DTRSM_D(side, uplo, transa, diag, m, n, alpha, alphad, +C + a, ad, lda, b, bd, ldb) +C IMPLICIT NONE +CC +CC -- Reference BLAS level3 routine -- +CC -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +CC -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +CC +CC .. Scalar Arguments .. +C DOUBLE PRECISION alpha +C DOUBLE PRECISION alphad +C INTEGER lda, ldb, m, n +C CHARACTER diag, side, transa, uplo +CC .. +CC .. Array Arguments .. +C DOUBLE PRECISION a(lda, *), b(ldb, *) +C DOUBLE PRECISION ad(lda, *), bd(ldb, *) +C EXTERNAL LSAME +CC .. +CC +CC ===================================================================== +CC +CC .. External Functions .. +C LOGICAL LSAME +CC .. +CC .. External Subroutines .. +C EXTERNAL XERBLA +CC .. +CC .. Intrinsic Functions .. +C INTRINSIC MAX +CC .. +CC .. Local Scalars .. +C DOUBLE PRECISION temp +C DOUBLE PRECISION tempd +C INTEGER i, info, j, k, nrowa +C LOGICAL lside, nounit, upper +CC .. +CC .. Parameters .. +C DOUBLE PRECISION zero +C PARAMETER (zero=0.0d+0) +C INTEGER max1 +C INTEGER max2 +C DOUBLE PRECISION temp0 +CC .. +CC +CC Test the input parameters. +CC +C lside = LSAME(side, 'L') +C IF (lside) THEN +C nrowa = m +C ELSE +C nrowa = n +C END IF +C nounit = LSAME(diag, 'N') +C upper = LSAME(uplo, 'U') +CC +C info = 0 +C IF (.NOT.lside .AND. (.NOT.LSAME(side, 'R'))) THEN +C info = 1 +C ELSE IF (.NOT.upper .AND. (.NOT.LSAME(uplo, 'L'))) THEN +C info = 2 +C ELSE IF (.NOT.LSAME(transa, 'N') .AND. (.NOT.LSAME(transa, 'T')) +C + .AND. (.NOT.LSAME(transa, 'C'))) THEN +C info = 3 +C ELSE IF (.NOT.LSAME(diag, 'U') .AND. (.NOT.LSAME(diag, 'N'))) THEN +C info = 4 +C ELSE IF (m .LT. 0) THEN +C info = 5 +C ELSE IF (n .LT. 0) THEN +C info = 6 +C ELSE +C IF (1 .LT. nrowa) THEN +C max1 = nrowa +C ELSE +C max1 = 1 +C END IF +C IF (lda .LT. max1) THEN +C info = 9 +C ELSE +C IF (1 .LT. m) THEN +C max2 = m +C ELSE +C max2 = 1 +C END IF +C IF (ldb .LT. max2) info = 11 +C END IF +C END IF +C IF (info .NE. 0) THEN +C CALL XERBLA('DTRSM ', info) +C RETURN +C ELSE IF (m .EQ. 0 .OR. n .EQ. 0) THEN +CC +CC Quick return if possible. +CC +C RETURN +C ELSE IF (alpha .EQ. zero) THEN +CC +CC And when alpha.eq.zero. +CC +C DO j=1,n +C DO i=1,m +C bd(i, j) = 0.D0 +C b(i, j) = zero +C ENDDO +C ENDDO +C RETURN +C ELSE +CC +CC Start the operations. +CC +C IF (lside) THEN +C IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*inv( A )*B. +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C bd(i, j) = b(i, j)*alphad + alpha*bd(i, j) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=m,1,-1 +C IF (nounit) THEN +C temp0 = b(k, j)/a(k, k) +C bd(k, j) = (bd(k, j)-temp0*ad(k, k))/a(k, k) +C b(k, j) = temp0 +C END IF +C DO i=1,k-1 +C bd(i, j) = bd(i, j) - a(i, k)*bd(k, j) - b(k, j)*ad( +C + i, k) +C b(i, j) = b(i, j) - b(k, j)*a(i, k) +C ENDDO +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=1,m +C bd(i, j) = b(i, j)*alphad + alpha*bd(i, j) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,m +C IF (nounit) THEN +C temp0 = b(k, j)/a(k, k) +C bd(k, j) = (bd(k, j)-temp0*ad(k, k))/a(k, k) +C b(k, j) = temp0 +C END IF +C DO i=k+1,m +C bd(i, j) = bd(i, j) - a(i, k)*bd(k, j) - b(k, j)*ad( +C + i, k) +C b(i, j) = b(i, j) - b(k, j)*a(i, k) +C ENDDO +C ENDDO +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*inv( A**T )*B. +CC +C DO j=1,n +C DO i=1,m +C tempd = b(i, j)*alphad + alpha*bd(i, j) +C temp = alpha*b(i, j) +C DO k=1,i-1 +C tempd = tempd - b(k, j)*ad(k, i) - a(k, i)*bd(k, j) +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C IF (nounit) THEN +C temp0 = temp/a(i, i) +C tempd = (tempd-temp0*ad(i, i))/a(i, i) +C temp = temp0 +C END IF +C bd(i, j) = tempd +C b(i, j) = temp +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=m,1,-1 +C tempd = b(i, j)*alphad + alpha*bd(i, j) +C temp = alpha*b(i, j) +C DO k=i+1,m +C tempd = tempd - b(k, j)*ad(k, i) - a(k, i)*bd(k, j) +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C IF (nounit) THEN +C temp0 = temp/a(i, i) +C tempd = (tempd-temp0*ad(i, i))/a(i, i) +C temp = temp0 +C END IF +C bd(i, j) = tempd +C b(i, j) = temp +C ENDDO +C ENDDO +C END IF +C ELSE IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*B*inv( A ). +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C bd(i, j) = b(i, j)*alphad + alpha*bd(i, j) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,j-1 +C DO i=1,m +C bd(i, j) = bd(i, j) - b(i, k)*ad(k, j) - a(k, j)*bd(i +C + , k) +C b(i, j) = b(i, j) - a(k, j)*b(i, k) +C ENDDO +C ENDDO +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, j)/a(j, j) +C bd(i, j) = (bd(i, j)-temp0*ad(j, j))/a(j, j) +C b(i, j) = temp0 +C ENDDO +C END IF +C ENDDO +C ELSE +C DO j=n,1,-1 +C DO i=1,m +C bd(i, j) = b(i, j)*alphad + alpha*bd(i, j) +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=j+1,n +C DO i=1,m +C bd(i, j) = bd(i, j) - b(i, k)*ad(k, j) - a(k, j)*bd(i +C + , k) +C b(i, j) = b(i, j) - a(k, j)*b(i, k) +C ENDDO +C ENDDO +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, j)/a(j, j) +C bd(i, j) = (bd(i, j)-temp0*ad(j, j))/a(j, j) +C b(i, j) = temp0 +C ENDDO +C END IF +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*B*inv( A**T ). +CC +C DO k=n,1,-1 +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, k)/a(k, k) +C bd(i, k) = (bd(i, k)-temp0*ad(k, k))/a(k, k) +C b(i, k) = temp0 +C ENDDO +C END IF +C DO j=1,k-1 +C DO i=1,m +C bd(i, j) = bd(i, j) - b(i, k)*ad(j, k) - a(j, k)*bd(i, k +C + ) +C b(i, j) = b(i, j) - a(j, k)*b(i, k) +C ENDDO +C ENDDO +C DO i=1,m +C bd(i, k) = b(i, k)*alphad + alpha*bd(i, k) +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C ELSE +C DO k=1,n +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, k)/a(k, k) +C bd(i, k) = (bd(i, k)-temp0*ad(k, k))/a(k, k) +C b(i, k) = temp0 +C ENDDO +C END IF +C DO j=k+1,n +C DO i=1,m +C bd(i, j) = bd(i, j) - b(i, k)*ad(j, k) - a(j, k)*bd(i, k +C + ) +C b(i, j) = b(i, j) - a(j, k)*b(i, k) +C ENDDO +C ENDDO +C DO i=1,m +C bd(i, k) = b(i, k)*alphad + alpha*bd(i, k) +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C END IF +CC +C RETURN +CC +CC End of DTRSM +CC +C END IF +C END SUBROUTINE DTRSM_D + + SUBROUTINE DTRSM_D(SIDE, UPLO, TRANSA, DIAG, M, N, + + ALPHA, ALPHAD, A, AD, LDA, B, BD, LDB) +C +C Forward-mode derivative of DTRSM via the black-box / Giles-style +C approach: uses ONLY the original, undifferentiated +C DTRSM and DTRMM -- no hand-differentiated loops. +C + IMPLICIT NONE + CHARACTER SIDE, UPLO, TRANSA, DIAG + INTEGER M, N, LDA, LDB + DOUBLE PRECISION ALPHA, ALPHAD + DOUBLE PRECISION A(LDA,*), AD(LDA,*) + DOUBLE PRECISION B(LDB,*), BD(LDB,*) + + DOUBLE PRECISION X(M,N), T(M,N) + INTEGER I, J + + IF (M.EQ.0 .OR. N.EQ.0) RETURN + +C Step 1: X = solution of the original forward call +C (private copy -- B itself must stay untouched) + DO J = 1, N + DO I = 1, M + X(I,J) = B(I,J) + END DO + END DO + CALL DTRSM(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, X, M) + +C Step 2: T = op(Ad)*X (SIDE='L') or X*op(Ad) (SIDE='R') + DO J = 1, N + DO I = 1, M + T(I,J) = X(I,J) + END DO + END DO + CALL DTRMM(SIDE, UPLO, TRANSA, DIAG, M, N, 1.0D0, AD, LDA, T, M) + +C Step 3: build RHS = alphad*B + alpha*Bd - T, in place in BD + DO J = 1, N + DO I = 1, M + BD(I,J) = ALPHAD*B(I,J) + ALPHA*BD(I,J) - T(I,J) + END DO + END DO + +C Step 4: solve the SAME triangular system again -> dX, overwrites BD + CALL DTRSM(SIDE, UPLO, TRANSA, DIAG, M, N, 1.0D0, A, LDA, BD, LDB) + + RETURN + END SUBROUTINE DTRSM_D + diff --git a/BLAS/src/dtrsm_dv.f b/BLAS/src/dtrsm_dv.f new file mode 100644 index 0000000..9c91a3e --- /dev/null +++ b/BLAS/src/dtrsm_dv.f @@ -0,0 +1,637 @@ +C Generated by TAPENADE (INRIA, Ecuador team) +C Tapenade 3.16 (develop) - 14 Mar 2025 11:54 +C +C Differentiation of dtrsm in forward (tangent) mode (with options multiDirectional): +C variations of useful results: b +C with respect to varying inputs: alpha a b +C RW status of diff variables: alpha:in a:in b:in-out +C> \brief \b DTRSM +C +C =========== DOCUMENTATION =========== +C +C Online html documentation available at +C http://www.netlib.org/lapack/explore-html/ +C +C Definition: +C =========== +C +C SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +C +C .. Scalar Arguments .. +C DOUBLE PRECISION ALPHA +C INTEGER LDA,LDB,M,N +C CHARACTER DIAG,SIDE,TRANSA,UPLO +C .. +C .. Array Arguments .. +C DOUBLE PRECISION A(LDA,*),B(LDB,*) +C .. +C +C +C> \par Purpose: +C ============= +C> +C> \verbatim +C> +C> DTRSM solves one of the matrix equations +C> +C> op( A )*X = alpha*B, or X*op( A ) = alpha*B, +C> +C> where alpha is a scalar, X and B are m by n matrices, A is a unit, or +C> non-unit, upper or lower triangular matrix and op( A ) is one of +C> +C> op( A ) = A or op( A ) = A**T. +C> +C> The matrix X is overwritten on B. +C> \endverbatim +C +C Arguments: +C ========== +C +C> \param[in] SIDE +C> \verbatim +C> SIDE is CHARACTER*1 +C> On entry, SIDE specifies whether op( A ) appears on the left +C> or right of X as follows: +C> +C> SIDE = 'L' or 'l' op( A )*X = alpha*B. +C> +C> SIDE = 'R' or 'r' X*op( A ) = alpha*B. +C> \endverbatim +C> +C> \param[in] UPLO +C> \verbatim +C> UPLO is CHARACTER*1 +C> On entry, UPLO specifies whether the matrix A is an upper or +C> lower triangular matrix as follows: +C> +C> UPLO = 'U' or 'u' A is an upper triangular matrix. +C> +C> UPLO = 'L' or 'l' A is a lower triangular matrix. +C> \endverbatim +C> +C> \param[in] TRANSA +C> \verbatim +C> TRANSA is CHARACTER*1 +C> On entry, TRANSA specifies the form of op( A ) to be used in +C> the matrix multiplication as follows: +C> +C> TRANSA = 'N' or 'n' op( A ) = A. +C> +C> TRANSA = 'T' or 't' op( A ) = A**T. +C> +C> TRANSA = 'C' or 'c' op( A ) = A**T. +C> \endverbatim +C> +C> \param[in] DIAG +C> \verbatim +C> DIAG is CHARACTER*1 +C> On entry, DIAG specifies whether or not A is unit triangular +C> as follows: +C> +C> DIAG = 'U' or 'u' A is assumed to be unit triangular. +C> +C> DIAG = 'N' or 'n' A is not assumed to be unit +C> triangular. +C> \endverbatim +C> +C> \param[in] M +C> \verbatim +C> M is INTEGER +C> On entry, M specifies the number of rows of B. M must be at +C> least zero. +C> \endverbatim +C> +C> \param[in] N +C> \verbatim +C> N is INTEGER +C> On entry, N specifies the number of columns of B. N must be +C> at least zero. +C> \endverbatim +C> +C> \param[in] ALPHA +C> \verbatim +C> ALPHA is DOUBLE PRECISION. +C> On entry, ALPHA specifies the scalar alpha. When alpha is +C> zero then A is not referenced and B need not be set before +C> entry. +C> \endverbatim +C> +C> \param[in] A +C> \verbatim +C> A is DOUBLE PRECISION array, dimension ( LDA, k ), +C> where k is m when SIDE = 'L' or 'l' +C> and k is n when SIDE = 'R' or 'r'. +C> Before entry with UPLO = 'U' or 'u', the leading k by k +C> upper triangular part of the array A must contain the upper +C> triangular matrix and the strictly lower triangular part of +C> A is not referenced. +C> Before entry with UPLO = 'L' or 'l', the leading k by k +C> lower triangular part of the array A must contain the lower +C> triangular matrix and the strictly upper triangular part of +C> A is not referenced. +C> Note that when DIAG = 'U' or 'u', the diagonal elements of +C> A are not referenced either, but are assumed to be unity. +C> \endverbatim +C> +C> \param[in] LDA +C> \verbatim +C> LDA is INTEGER +C> On entry, LDA specifies the first dimension of A as declared +C> in the calling (sub) program. When SIDE = 'L' or 'l' then +C> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +C> then LDA must be at least max( 1, n ). +C> \endverbatim +C> +C> \param[in,out] B +C> \verbatim +C> B is DOUBLE PRECISION array, dimension ( LDB, N ) +C> Before entry, the leading m by n part of the array B must +C> contain the right-hand side matrix B, and on exit is +C> overwritten by the solution matrix X. +C> \endverbatim +C> +C> \param[in] LDB +C> \verbatim +C> LDB is INTEGER +C> On entry, LDB specifies the first dimension of B as declared +C> in the calling (sub) program. LDB must be at least +C> max( 1, m ). +C> \endverbatim +C +C Authors: +C ======== +C +C> \author Univ. of Tennessee +C> \author Univ. of California Berkeley +C> \author Univ. of Colorado Denver +C> \author NAG Ltd. +C +C> \ingroup trsm +C +C> \par Further Details: +C ===================== +C> +C> \verbatim +C> +C> Level 3 Blas routine. +C> +C> +C> -- Written on 8-February-1989. +C> Jack Dongarra, Argonne National Laboratory. +C> Iain Duff, AERE Harwell. +C> Jeremy Du Croz, Numerical Algorithms Group Ltd. +C> Sven Hammarling, Numerical Algorithms Group Ltd. +C> \endverbatim +C> +C ===================================================================== +C SUBROUTINE DTRSM_DV(side, uplo, transa, diag, m, n, alpha, alphad +C + , a, ad, lda, b, bd, ldb, nbdirs) +C IMPLICIT NONE +C INCLUDE 'DIFFSIZES.inc' +CC Hint: nbdirsmax should be the maximum number of differentiation directions +CC +CC -- Reference BLAS level3 routine -- +CC -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +CC -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +CC +CC .. Scalar Arguments .. +C DOUBLE PRECISION alpha +C DOUBLE PRECISION alphad(nbdirsmax) +C INTEGER lda, ldb, m, n +C CHARACTER diag, side, transa, uplo +CC .. +CC .. Array Arguments .. +C DOUBLE PRECISION a(lda, *), b(ldb, *) +C DOUBLE PRECISION ad(nbdirsmax, lda, *), bd(nbdirsmax, ldb, *) +C EXTERNAL LSAME +CC .. +CC +CC ===================================================================== +CC +CC .. External Functions .. +C LOGICAL LSAME +CC .. +CC .. External Subroutines .. +C EXTERNAL XERBLA +CC .. +CC .. Intrinsic Functions .. +C INTRINSIC MAX +CC .. +CC .. Local Scalars .. +C DOUBLE PRECISION temp +C DOUBLE PRECISION tempd(nbdirsmax) +C INTEGER i, info, j, k, nrowa +C LOGICAL lside, nounit, upper +CC .. +CC .. Parameters .. +C DOUBLE PRECISION zero +C PARAMETER (zero=0.0d+0) +C INTEGER max1 +C INTEGER max2 +C INTEGER nd +C DOUBLE PRECISION temp0 +C INTEGER nbdirs +CC .. +CC +CC Test the input parameters. +CC +CC Check 0 < nbdirs <= nbdirsmax (required by DIFFSIZES.inc) +C IF (nbdirs.LE.0 .OR. nbdirs.GT.nbdirsmax) THEN +C WRITE(*,'(A,I0,A,I0,A)') 'Error: nbdirs=', nbdirs, +C + ' must be in 1..nbdirsmax=', nbdirsmax, '. Stopping.' +C STOP 1 +C END IF +CC +C lside = LSAME(side, 'L') +C IF (lside) THEN +C nrowa = m +C ELSE +C nrowa = n +C END IF +C nounit = LSAME(diag, 'N') +C upper = LSAME(uplo, 'U') +CC +C info = 0 +C IF (.NOT.lside .AND. (.NOT.LSAME(side, 'R'))) THEN +C info = 1 +C ELSE IF (.NOT.upper .AND. (.NOT.LSAME(uplo, 'L'))) THEN +C info = 2 +C ELSE IF (.NOT.LSAME(transa, 'N') .AND. (.NOT.LSAME(transa, 'T')) +C + .AND. (.NOT.LSAME(transa, 'C'))) THEN +C info = 3 +C ELSE IF (.NOT.LSAME(diag, 'U') .AND. (.NOT.LSAME(diag, 'N'))) THEN +C info = 4 +C ELSE IF (m .LT. 0) THEN +C info = 5 +C ELSE IF (n .LT. 0) THEN +C info = 6 +C ELSE +C IF (1 .LT. nrowa) THEN +C max1 = nrowa +C ELSE +C max1 = 1 +C END IF +C IF (lda .LT. max1) THEN +C info = 9 +C ELSE +C IF (1 .LT. m) THEN +C max2 = m +C ELSE +C max2 = 1 +C END IF +C IF (ldb .LT. max2) info = 11 +C END IF +C END IF +C IF (info .NE. 0) THEN +C CALL XERBLA('DTRSM ', info) +C RETURN +C ELSE IF (m .EQ. 0 .OR. n .EQ. 0) THEN +CC +CC Quick return if possible. +CC +C RETURN +C ELSE IF (alpha .EQ. zero) THEN +CC +CC And when alpha.eq.zero. +CC +C DO j=1,n +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = 0.D0 +C ENDDO +C b(i, j) = zero +C ENDDO +C ENDDO +C RETURN +C ELSE +CC +CC Start the operations. +CC +C IF (lside) THEN +C IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*inv( A )*B. +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = b(i, j)*alphad(nd) + alpha*bd(nd, i, +C + j) +C ENDDO +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=m,1,-1 +C IF (nounit) THEN +C temp0 = b(k, j)/a(k, k) +C DO nd=1,nbdirs +C bd(nd, k, j) = (bd(nd, k, j)-temp0*ad(nd, k, k))/a +C + (k, k) +C ENDDO +C b(k, j) = temp0 +C END IF +C DO i=1,k-1 +C DO nd=1,nbdirs +C bd(nd, i, j) = bd(nd, i, j) - a(i, k)*bd(nd, k, j) +C + - b(k, j)*ad(nd, i, k) +C ENDDO +C b(i, j) = b(i, j) - b(k, j)*a(i, k) +C ENDDO +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = b(i, j)*alphad(nd) + alpha*bd(nd, i, +C + j) +C ENDDO +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,m +C IF (nounit) THEN +C temp0 = b(k, j)/a(k, k) +C DO nd=1,nbdirs +C bd(nd, k, j) = (bd(nd, k, j)-temp0*ad(nd, k, k))/a +C + (k, k) +C ENDDO +C b(k, j) = temp0 +C END IF +C DO i=k+1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = bd(nd, i, j) - a(i, k)*bd(nd, k, j) +C + - b(k, j)*ad(nd, i, k) +C ENDDO +C b(i, j) = b(i, j) - b(k, j)*a(i, k) +C ENDDO +C ENDDO +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*inv( A**T )*B. +CC +C DO j=1,n +C DO i=1,m +C DO nd=1,nbdirs +C tempd(nd) = b(i, j)*alphad(nd) + alpha*bd(nd, i, j) +C ENDDO +C temp = alpha*b(i, j) +C DO k=1,i-1 +C DO nd=1,nbdirs +C tempd(nd) = tempd(nd) - b(k, j)*ad(nd, k, i) - a(k, +C + i)*bd(nd, k, j) +C ENDDO +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C IF (nounit) THEN +C temp0 = temp/a(i, i) +C DO nd=1,nbdirs +C tempd(nd) = (tempd(nd)-temp0*ad(nd, i, i))/a(i, i) +C ENDDO +C temp = temp0 +C END IF +C DO nd=1,nbdirs +C bd(nd, i, j) = tempd(nd) +C ENDDO +C b(i, j) = temp +C ENDDO +C ENDDO +C ELSE +C DO j=1,n +C DO i=m,1,-1 +C DO nd=1,nbdirs +C tempd(nd) = b(i, j)*alphad(nd) + alpha*bd(nd, i, j) +C ENDDO +C temp = alpha*b(i, j) +C DO k=i+1,m +C DO nd=1,nbdirs +C tempd(nd) = tempd(nd) - b(k, j)*ad(nd, k, i) - a(k, +C + i)*bd(nd, k, j) +C ENDDO +C temp = temp - a(k, i)*b(k, j) +C ENDDO +C IF (nounit) THEN +C temp0 = temp/a(i, i) +C DO nd=1,nbdirs +C tempd(nd) = (tempd(nd)-temp0*ad(nd, i, i))/a(i, i) +C ENDDO +C temp = temp0 +C END IF +C DO nd=1,nbdirs +C bd(nd, i, j) = tempd(nd) +C ENDDO +C b(i, j) = temp +C ENDDO +C ENDDO +C END IF +C ELSE IF (LSAME(transa, 'N')) THEN +CC +CC Form B := alpha*B*inv( A ). +CC +C IF (upper) THEN +C DO j=1,n +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = b(i, j)*alphad(nd) + alpha*bd(nd, i, j) +C ENDDO +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=1,j-1 +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = bd(nd, i, j) - b(i, k)*ad(nd, k, j) - +C + a(k, j)*bd(nd, i, k) +C ENDDO +C b(i, j) = b(i, j) - a(k, j)*b(i, k) +C ENDDO +C ENDDO +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, j)/a(j, j) +C DO nd=1,nbdirs +C bd(nd, i, j) = (bd(nd, i, j)-temp0*ad(nd, j, j))/a(j +C + , j) +C ENDDO +C b(i, j) = temp0 +C ENDDO +C END IF +C ENDDO +C ELSE +C DO j=n,1,-1 +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = b(i, j)*alphad(nd) + alpha*bd(nd, i, j) +C ENDDO +C b(i, j) = alpha*b(i, j) +C ENDDO +C DO k=j+1,n +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = bd(nd, i, j) - b(i, k)*ad(nd, k, j) - +C + a(k, j)*bd(nd, i, k) +C ENDDO +C b(i, j) = b(i, j) - a(k, j)*b(i, k) +C ENDDO +C ENDDO +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, j)/a(j, j) +C DO nd=1,nbdirs +C bd(nd, i, j) = (bd(nd, i, j)-temp0*ad(nd, j, j))/a(j +C + , j) +C ENDDO +C b(i, j) = temp0 +C ENDDO +C END IF +C ENDDO +C END IF +C ELSE IF (upper) THEN +CC +CC Form B := alpha*B*inv( A**T ). +CC +C DO k=n,1,-1 +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, k)/a(k, k) +C DO nd=1,nbdirs +C bd(nd, i, k) = (bd(nd, i, k)-temp0*ad(nd, k, k))/a(k, +C + k) +C ENDDO +C b(i, k) = temp0 +C ENDDO +C END IF +C DO j=1,k-1 +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = bd(nd, i, j) - b(i, k)*ad(nd, j, k) - a +C + (j, k)*bd(nd, i, k) +C ENDDO +C b(i, j) = b(i, j) - a(j, k)*b(i, k) +C ENDDO +C ENDDO +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, k) = b(i, k)*alphad(nd) + alpha*bd(nd, i, k) +C ENDDO +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C ELSE +C DO k=1,n +C IF (nounit) THEN +C DO i=1,m +C temp0 = b(i, k)/a(k, k) +C DO nd=1,nbdirs +C bd(nd, i, k) = (bd(nd, i, k)-temp0*ad(nd, k, k))/a(k, +C + k) +C ENDDO +C b(i, k) = temp0 +C ENDDO +C END IF +C DO j=k+1,n +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, j) = bd(nd, i, j) - b(i, k)*ad(nd, j, k) - a +C + (j, k)*bd(nd, i, k) +C ENDDO +C b(i, j) = b(i, j) - a(j, k)*b(i, k) +C ENDDO +C ENDDO +C DO i=1,m +C DO nd=1,nbdirs +C bd(nd, i, k) = b(i, k)*alphad(nd) + alpha*bd(nd, i, k) +C ENDDO +C b(i, k) = alpha*b(i, k) +C ENDDO +C ENDDO +C END IF +CC +C RETURN +CC +CC End of DTRSM +CC +C END IF +C END + + SUBROUTINE DTRSM_DV(SIDE, UPLO, TRANSA, DIAG, M, N, + + ALPHA, ALPHAD, A, AD, LDA, B, BD, LDB, NBDIRS) +C +C Vector forward-mode derivative of DTRSM (black-box / Giles-style): +C differentiates op(A)*X = alpha*B directly, using only the original +C DTRSM and DTRMM. Self-contained, does not call DTRSM_D. +C +C AD/BD store direction as the FASTEST-varying index (nbdirsmax, LDA/LDB, *), +C so each direction's slice is gathered into a plain work array before use. +C + IMPLICIT NONE + INCLUDE 'DIFFSIZES.inc' + CHARACTER SIDE, UPLO, TRANSA, DIAG + INTEGER M, N, LDA, LDB, NBDIRS + DOUBLE PRECISION ALPHA + DOUBLE PRECISION ALPHAD(NBDIRSMAX) + DOUBLE PRECISION A(LDA,*), AD(NBDIRSMAX,LDA,*) + DOUBLE PRECISION B(LDB,*), BD(NBDIRSMAX,LDB,*) + + DOUBLE PRECISION X(M,N), T(M,N) + DOUBLE PRECISION AD_SLICE(LDA,N), BD_SLICE(LDB,N) + INTEGER ND, I, J + + IF (NBDIRS.LE.0 .OR. NBDIRS.GT.NBDIRSMAX) THEN + WRITE(*,'(A,I0,A,I0,A)') 'Error: nbdirs=', NBDIRS, + + ' must be in 1..nbdirsmax=', NBDIRSMAX, '. Stopping.' + STOP 1 + END IF + + IF (M.EQ.0 .OR. N.EQ.0) RETURN + +C X = solution of the original forward call (same for every direction) + DO J = 1, N + DO I = 1, M + X(I,J) = B(I,J) + END DO + END DO + CALL DTRSM(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, X, M) + + DO ND = 1, NBDIRS + +C Gather this direction out of the interleaved AD/BD layout + DO J = 1, N + DO I = 1, M + AD_SLICE(I,J) = AD(ND,I,J) + BD_SLICE(I,J) = BD(ND,I,J) + END DO + END DO + +C T = op(Ad)*X (SIDE='L') or X*op(Ad) (SIDE='R') + DO J = 1, N + DO I = 1, M + T(I,J) = X(I,J) + END DO + END DO + CALL DTRMM(SIDE, UPLO, TRANSA, DIAG, M, N, 1.0D0, + + AD_SLICE, LDA, T, M) + +C RHS = alphad*B + alpha*Bd - T, built in place in BD_SLICE + DO J = 1, N + DO I = 1, M + BD_SLICE(I,J) = ALPHAD(ND)*B(I,J) + ALPHA*BD_SLICE(I,J) + + - T(I,J) + END DO + END DO + +C Solve the same triangular system again -> dX for this direction + CALL DTRSM(SIDE, UPLO, TRANSA, DIAG, M, N, 1.0D0, + + A, LDA, BD_SLICE, LDB) + +C Scatter the result back into the interleaved layout + DO J = 1, N + DO I = 1, M + BD(ND,I,J) = BD_SLICE(I,J) + END DO + END DO + + END DO + + RETURN + END SUBROUTINE DTRSM_DV diff --git a/BLAS/test/test_dtrsm.f90 b/BLAS/test/test_dtrsm.f90 new file mode 100644 index 0000000..11ada7b --- /dev/null +++ b/BLAS/test/test_dtrsm.f90 @@ -0,0 +1,201 @@ +! Test program for DTRSM differentiation +! Generated automatically by run_tapenade_blas.py +! Using REAL*8 precision + +program test_dtrsm + implicit none + + integer :: seed_array(33) + + external :: dtrsm + external :: dtrsm_d + + ! Test parameters + integer, parameter :: n = 4 ! Matrix/vector size for test + integer, parameter :: max_size = n ! Maximum array dimension (rows/cols of matrices) + integer, parameter :: lda = max_size, ldb = max_size, ldc = max_size ! Leading dimensions + + character :: side + character :: uplo + character :: transa + character :: diag + integer :: msize + integer :: nsize + real(8) :: alpha + real(8), dimension(max_size,max_size) :: a + integer :: lda_val + real(8), dimension(max_size,max_size) :: b + integer :: ldb_val + + ! Derivative variables + real(8) :: alpha_d + real(8), dimension(max_size,max_size) :: a_d + real(8), dimension(max_size,max_size) :: b_d + + ! Storage variables for inout parameters + real(8), dimension(max_size,max_size) :: b_output + + ! Array restoration variables for numerical differentiation + real(8) :: alpha_orig + real(8), dimension(max_size,max_size) :: b_orig + real(8), dimension(max_size,max_size) :: a_orig + + ! Variables for central difference computation + real(8), dimension(max_size,max_size) :: b_forward, b_backward + ! Scalar variables for central difference computation + real(8) :: central_diff, ad_result + logical :: has_large_errors + + ! Variables for storing original derivative values + real(8) :: alpha_d_orig + real(8), dimension(max_size,max_size) :: b_d_orig + real(8), dimension(max_size,max_size) :: a_d_orig + + ! Temporary variables for matrix initialization + real(4) :: temp_real, temp_imag + integer :: i, j + + ! Initialize random seed for reproducible results + seed_array = 42 + call random_seed(put=seed_array) + + side = 'L' + uplo = 'U' + transa = 'N' + diag = 'N' + msize = n + nsize = n + call random_number(alpha) + alpha = alpha * 2.0d0 - 1.0d0 ! Scale to [-1,1] + call random_number(a) + a = a * 2.0d0 - 1.0d0 ! Scale to [-1,1] + lda_val = lda ! LDA must be at least max( 1 + call random_number(b) + b = b * 2.0d0 - 1.0d0 ! Scale to [-1,1] + ldb_val = ldb + + ! Initialize input derivatives to random values + call random_number(alpha_d) + alpha_d = alpha_d * 2.0e0 - 1.0e0 ! Scale to [-1,1] + call random_number(b_d) + b_d = b_d * 2.0e0 - 1.0e0 ! Scale to [-1,1] + call random_number(a_d) + a_d = a_d * 2.0e0 - 1.0e0 ! Scale to [-1,1] + + ! Store initial derivative values after random initialization + alpha_d_orig = alpha_d + b_d_orig = b_d + a_d_orig = a_d + + ! Store original values for central difference computation + alpha_orig = alpha + b_orig = b + a_orig = a + + write(*,*) 'Testing DTRSM' + ! Store input values of inout parameters before first function call + b_orig = b + + ! Re-initialize data for differentiated function + ! Only reinitialize inout parameters - keep input-only parameters unchanged + + ! side already has correct value from original call + ! uplo already has correct value from original call + ! transa already has correct value from original call + ! diag already has correct value from original call + msize = n + nsize = n + ! alpha already has correct value from original call + ! a already has correct value from original call + lda_val = lda ! LDA must be at least max( 1 + b = b_orig + ldb_val = ldb + + ! Call the differentiated function + call dtrsm_d(side, uplo, transa, diag, msize, nsize, alpha, alpha_d, a, a_d, lda_val, b, b_d, ldb_val) + + ! Print results and compare + write(*,*) 'Function calls completed successfully' + + ! Numerical differentiation check + call check_derivatives_numerically() + + write(*,*) 'Test completed successfully' + +contains + + subroutine check_derivatives_numerically() + implicit none + real(8), parameter :: h = 1.0e-6 ! Step size for finite differences + real(8) :: relative_error, max_error + real(8) :: output_orig, output_pert + real(8) :: numerical_result, analytical_result + real(8) :: abs_error, abs_reference, error_bound + integer :: i, j + + max_error = 0.0e0 + has_large_errors = .false. + + write(*,*) 'Checking derivatives against numerical differentiation:' + write(*,*) 'Step size h =', h + + ! Tolerance thresholds: rtol=1.0e-5, atol=1.0e-5 + + ! Original values already stored in main program + + ! Central difference computation: f(x + h) - f(x - h) / (2h) + ! Forward perturbation: f(x + h) + alpha = alpha_orig + h * alpha_d_orig + b = b_orig + h * b_d_orig + a = a_orig + h * a_d_orig + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + ! Store forward perturbation results + b_forward = b + + ! Backward perturbation: f(x - h) + alpha = alpha_orig - h * alpha_d_orig + b = b_orig - h * b_d_orig + a = a_orig - h * a_d_orig + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + ! Store backward perturbation results + b_backward = b + + ! Compute central differences and compare with AD results + ! Check derivatives for output B + do j = 1, min(2, n) ! Check only first few elements + do i = 1, min(2, n) + ! Central difference: (f(x+h) - f(x-h)) / (2h) + central_diff = (b_forward(i,j) - b_backward(i,j)) / (2.0e0 * h) + ! AD result + ad_result = b_d(i,j) + ! Error check: |a - b| > atol + rtol * |b| + abs_error = abs(central_diff - ad_result) + abs_reference = abs(ad_result) + error_bound = 1.0e-5 + 1.0e-5 * abs_reference + if (abs_error > error_bound) then + has_large_errors = .true. + relative_error = abs_error / max(abs_reference, 1.0e-10) + write(*,*) 'Large error in output B(', i, ',', j, '):' + write(*,*) ' Central diff: ', central_diff + write(*,*) ' AD result: ', ad_result + write(*,*) ' Absolute error:', abs_error + write(*,*) ' Error bound:', error_bound + write(*,*) ' Relative error:', relative_error + end if + ! Track max error for reporting (normalized) + relative_error = abs_error / max(abs_reference, 1.0e-10) + max_error = max(max_error, relative_error) + end do + end do + + write(*,*) 'Maximum relative error:', max_error + write(*,*) 'Tolerance thresholds: rtol=1.0e-5, atol=1.0e-5' + if (has_large_errors) then + write(*,*) 'FAIL: Derivatives are outside tolerance' + else + write(*,*) 'PASS: Derivatives are within tolerance (rtol + atol)' + end if + + end subroutine check_derivatives_numerically + +end program test_dtrsm diff --git a/BLAS/test/test_dtrsm_reverse.f90 b/BLAS/test/test_dtrsm_reverse.f90 new file mode 100644 index 0000000..ebe5cb7 --- /dev/null +++ b/BLAS/test/test_dtrsm_reverse.f90 @@ -0,0 +1,254 @@ +! Test program for DTRSM reverse mode (adjoint) differentiation +! Generated automatically by run_tapenade_blas.py +! Using REAL*8 precision +! Verification uses VJP methodology with finite differences + +program test_dtrsm_reverse + implicit none + + integer :: seed_array(33) + + external :: dtrsm + external :: dtrsm_b + + ! Test parameters + integer, parameter :: n = 4 ! Matrix/vector size for test + integer, parameter :: max_size = n ! Maximum array dimension (rows/cols of matrices) + integer, parameter :: lda = max_size, ldb = max_size, ldc = max_size ! Leading dimensions + + character :: side + character :: uplo + character :: transa + character :: diag + integer :: msize + integer :: nsize + real(8) :: alpha + real(8), dimension(max_size,max_size) :: a + integer :: lda_val + real(8), dimension(max_size,max_size) :: b + integer :: ldb_val + + ! Adjoint variables (reverse mode) + ! In reverse mode: output adjoints are INPUT (cotangents/seeds) + ! input adjoints are OUTPUT (computed gradients) + real(8) :: alphab + real(8), dimension(max_size,max_size) :: ab + real(8), dimension(max_size,max_size) :: bb + + ! Storage for original values (for VJP verification) + real(8) :: alpha_orig + real(8), dimension(max_size,max_size) :: a_orig + real(8), dimension(max_size,max_size) :: b_orig + + ! Variables for VJP verification via finite differences + real(8), dimension(max_size,max_size) :: b_plus, b_minus + + ! Saved cotangents (output adjoints) for VJP verification + real(8), dimension(max_size,max_size) :: bb_orig + real(8), parameter :: h = 1.0e-7 + real(8) :: vjp_ad, vjp_fd, relative_error, max_error, abs_error, abs_reference, error_bound + logical :: has_large_errors + integer :: i, j + real(8), dimension(max_size*max_size) :: temp_products ! For sorted summation + integer :: n_products + + ! Initialize random seed for reproducibility + seed_array = 42 + call random_seed(put=seed_array) + + ! Initialize primal values + side = 'L' + uplo = 'U' + transa = 'N' + diag = 'N' + msize = n + nsize = n + call random_number(alpha) + alpha = alpha * 2.0d0 - 1.0d0 + call random_number(a) + a = a * 2.0d0 - 1.0d0 + lda_val = lda + call random_number(b) + b = b * 2.0d0 - 1.0d0 + ldb_val = ldb + + ! Store original primal values + alpha_orig = alpha + a_orig = a + b_orig = b + + write(*,*) 'Testing DTRSM' + + ! Initialize output adjoints (cotangents) with random values + ! These are the 'seeds' for reverse mode + call random_number(bb) + bb = bb * 2.0d0 - 1.0d0 + + ! Save output adjoints (cotangents) for VJP verification + ! Note: output adjoints may be modified by reverse mode function + bb_orig = bb + + ! Initialize input adjoints to zero (they will be computed) + ab = 0.0d0 + alphab = 0.0d0 + + ! Set ISIZE globals required by differentiated routine (dimension 2 of arrays). + ! Differentiated code checks they are set via check_ISIZE*_initialized. + call set_ISIZE2OFA(max_size) + + ! Call reverse mode differentiated function + call dtrsm_b(side, uplo, transa, diag, msize, nsize, alpha, alphab, a, ab, lda_val, b, bb, ldb_val) + + ! Reset ISIZE globals to uninitialized (-1) for completeness + call set_ISIZE2OFA(-1) + + ! VJP Verification using finite differences + ! For reverse mode, we verify: cotangent^T @ J @ direction = direction^T @ adjoint + ! Equivalently: cotangent^T @ (f(x+h*dir) - f(x-h*dir))/(2h) should equal dir^T @ computed_adjoint + call check_vjp_numerically() + + write(*,*) '' + write(*,*) 'Test completed successfully' + +contains + + subroutine check_vjp_numerically() + implicit none + + ! Direction vectors for VJP testing (like tangents in forward mode) + real(8) :: alpha_dir + real(8), dimension(max_size,max_size) :: a_dir + real(8), dimension(max_size,max_size) :: b_dir + + real(8), dimension(max_size,max_size) :: b_central_diff + + max_error = 0.0d0 + has_large_errors = .false. + + write(*,*) 'Function calls completed successfully' + + write(*,*) 'Checking derivatives against numerical differentiation:' + write(*,*) 'Step size h =', h + + ! Initialize random direction vectors for all inputs + call random_number(alpha_dir) + alpha_dir = alpha_dir * 2.0d0 - 1.0d0 + call random_number(a_dir) + a_dir = a_dir * 2.0d0 - 1.0d0 + call random_number(b_dir) + b_dir = b_dir * 2.0d0 - 1.0d0 + + ! Forward perturbation: f(x + h*dir) + alpha = alpha_orig + h * alpha_dir + a = a_orig + h * a_dir + b = b_orig + h * b_dir + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + b_plus = b + + ! Backward perturbation: f(x - h*dir) + alpha = alpha_orig - h * alpha_dir + a = a_orig - h * a_dir + b = b_orig - h * b_dir + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + b_minus = b + + ! Compute central differences: (f(x+h*dir) - f(x-h*dir)) / (2h) + b_central_diff = (b_plus - b_minus) / (2.0d0 * h) + + ! VJP verification: + ! cotangent^T @ central_diff should equal direction^T @ computed_adjoint + ! Left side: cotangent^T @ Jacobian @ direction (via finite differences, with sorted summation) + vjp_fd = 0.0d0 + ! Compute and sort products for b (FD) + n_products = 0 + do j = 1, n + do i = 1, n + n_products = n_products + 1 + temp_products(n_products) = bb_orig(i,j) * b_central_diff(i,j) + end do + end do + call sort_array(temp_products, n_products) + do i = 1, n_products + vjp_fd = vjp_fd + temp_products(i) + end do + + ! Right side: direction^T @ computed_adjoint (with sorted summation) + ! For INOUT parameters: use cb directly (it contains the computed input adjoint after reverse pass) + ! For pure inputs: use adjoint directly + vjp_ad = 0.0d0 + vjp_ad = vjp_ad + alpha_dir * alphab + ! Compute and sort products for a + n_products = 0 + do j = 1, n + do i = 1, n + n_products = n_products + 1 + temp_products(n_products) = a_dir(i,j) * ab(i,j) + end do + end do + call sort_array(temp_products, n_products) + do i = 1, n_products + vjp_ad = vjp_ad + temp_products(i) + end do + ! Compute and sort products for b + n_products = 0 + do j = 1, n + do i = 1, n + n_products = n_products + 1 + temp_products(n_products) = b_dir(i,j) * bb(i,j) + end do + end do + call sort_array(temp_products, n_products) + do i = 1, n_products + vjp_ad = vjp_ad + temp_products(i) + end do + + ! Error check: |vjp_fd - vjp_ad| > atol + rtol * |vjp_ad| + abs_error = abs(vjp_fd - vjp_ad) + abs_reference = abs(vjp_ad) + error_bound = 1.0e-5 + 1.0e-5 * abs_reference + if (abs_error > error_bound) then + has_large_errors = .true. + end if + + + if (abs_reference > 1.0e-10) then + relative_error = abs_error / abs_reference + else + relative_error = abs_error + end if + max_error = relative_error + + write(*,*) 'Maximum relative error:', max_error + write(*,*) 'Tolerance thresholds: rtol=1.0e-5, atol=1.0e-5' + if (has_large_errors) then + write(*,*) 'FAIL: Derivatives are outside tolerance' + else + write(*,*) 'PASS: Derivatives are within tolerance (rtol + atol)' + end if + + end subroutine check_vjp_numerically + + subroutine sort_array(arr, n) + implicit none + integer, intent(in) :: n + real(8), dimension(n), intent(inout) :: arr + integer :: i, j, min_idx + real(8) :: temp + + ! Simple selection sort + do i = 1, n-1 + min_idx = i + do j = i+1, n + if (abs(arr(j)) < abs(arr(min_idx))) then + min_idx = j + end if + end do + if (min_idx /= i) then + temp = arr(i) + arr(i) = arr(min_idx) + arr(min_idx) = temp + end if + end do + end subroutine sort_array + +end program test_dtrsm_reverse diff --git a/BLAS/test/test_dtrsm_vector_forward.f90 b/BLAS/test/test_dtrsm_vector_forward.f90 new file mode 100644 index 0000000..5e8404a --- /dev/null +++ b/BLAS/test/test_dtrsm_vector_forward.f90 @@ -0,0 +1,180 @@ +! Test program for DTRSM vector forward mode differentiation +! Generated automatically by run_tapenade_blas.py +! Using REAL*8 precision with nbdirsmax=4 + +program test_dtrsm_vector_forward + implicit none + include 'DIFFSIZES.inc' + + external :: dtrsm + external :: dtrsm_dv + + ! Test parameters + integer, parameter :: n = 4 ! Matrix/vector size for test + integer, parameter :: max_size = n ! Maximum array dimension + integer, parameter :: lda = max_size, ldb = max_size, ldc = max_size ! Leading dimensions + integer :: i, j, idir ! Loop counters + integer :: seed_array(33) ! Random seed + real(4) :: temp_real, temp_imag ! Temporary variables for complex initialization + + character :: side + character :: uplo + character :: transa + character :: diag + integer :: msize + integer :: nsize + real(8) :: alpha + real(8), dimension(max_size,max_size) :: a + integer :: lda_val + real(8), dimension(max_size,max_size) :: b + integer :: ldb_val + + ! Vector mode derivative variables (type-promoted) + ! Scalars become arrays(nbdirsmax), arrays gain extra dimension + real(8), dimension(nbdirsmax) :: alpha_dv + real(8), dimension(nbdirsmax,max_size,max_size) :: a_dv + real(8), dimension(nbdirsmax,max_size,max_size) :: b_dv + ! Declare variables for storing original values + real(8) :: alpha_orig + real(8), dimension(nbdirsmax) :: alpha_dv_orig + real(8), dimension(max_size,max_size) :: a_orig + real(8), dimension(nbdirsmax,max_size,max_size) :: a_dv_orig + real(8), dimension(max_size,max_size) :: b_orig + real(8), dimension(nbdirsmax,max_size,max_size) :: b_dv_orig + + seed_array = 42 + call random_seed(put=seed_array) + + + ! Initialize test parameters + msize = n + nsize = n + lda_val = lda + ldb_val = ldb + + ! Initialize test data with random numbers + ! Initialize random seed for reproducible results + seed_array = 42 + call random_seed(put=seed_array) + + side = 'L' + uplo = 'U' + transa = 'N' + diag = 'N' + call random_number(alpha) + alpha = alpha * 2.0d0 - 1.0d0 ! Scale to [-1,1] + call random_number(a) + a = a * 2.0d0 - 1.0d0 ! Scale to [-1,1] + call random_number(b) + b = b * 2.0d0 - 1.0d0 ! Scale to [-1,1] + + ! Initialize input derivatives to random values (exactly like scalar mode) + do idir = 1, nbdirsmax + call random_number(temp_real) + alpha_dv(idir) = temp_real * 2.0d0 - 1.0d0 + end do + do idir = 1, nbdirsmax + call random_number(a_dv(idir,:,:)) + a_dv(idir,:,:) = a_dv(idir,:,:) * 2.0d0 - 1.0d0 + end do + do idir = 1, nbdirsmax + call random_number(b_dv(idir,:,:)) + b_dv(idir,:,:) = b_dv(idir,:,:) * 2.0d0 - 1.0d0 + end do + + write(*,*) 'Testing DTRSM (Vector Forward Mode)' + ! Store original values before any function calls (critical for INOUT parameters) + alpha_orig = alpha + alpha_dv_orig = alpha_dv + a_orig = a + a_dv_orig = a_dv + b_orig = b + b_dv_orig = b_dv + + ! Call the vector mode differentiated function + + call dtrsm_dv(side, uplo, transa, diag, msize, nsize, alpha, alpha_dv, a, a_dv, lda_val, b, b_dv, ldb_val, nbdirsmax) + + ! Print results and compare + write(*,*) 'Function calls completed successfully' + + ! Numerical differentiation check + call check_derivatives_numerically() + + write(*,*) 'Vector forward mode test completed successfully' + +contains + + subroutine check_derivatives_numerically() + implicit none + real(8), parameter :: h = 1.0e-7 ! Step size for finite differences + real(8) :: relative_error, max_error + real(8) :: abs_error, abs_reference, error_bound + real(8) :: central_diff, ad_result + integer :: i, j, idir + logical :: has_large_errors + real(8), dimension(max_size,max_size) :: b_forward, b_backward + + max_error = 0.0e0 + has_large_errors = .false. + + write(*,*) 'Checking derivatives against numerical differentiation:' + write(*,*) 'Step size h =', h + write(*,*) 'Number of directions:', nbdirsmax + + ! Test each derivative direction separately + do idir = 1, nbdirsmax + + ! Forward perturbation: f(x + h * direction) + alpha = alpha_orig + h * alpha_dv_orig(idir) + a = a_orig + h * a_dv_orig(idir,:,:) + b = b_orig + h * b_dv_orig(idir,:,:) + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + b_forward = b + + ! Backward perturbation: f(x - h * direction) + alpha = alpha_orig - h * alpha_dv_orig(idir) + a = a_orig - h * a_dv_orig(idir,:,:) + b = b_orig - h * b_dv_orig(idir,:,:) + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + b_backward = b + + ! Compute central differences and compare with AD results + do j = 1, min(2, nsize) ! Check only first few elements + do i = 1, min(2, nsize) + ! Central difference: (f(x+h) - f(x-h)) / (2h) + central_diff = (b_forward(i,j) - b_backward(i,j)) / (2.0e0 * h) + ! AD result + ad_result = b_dv(idir,i,j) + ! Error check: |a - b| > atol + rtol * |b| + abs_error = abs(central_diff - ad_result) + abs_reference = abs(ad_result) + error_bound = 1.0e-5 + 1.0e-5 * abs_reference + if (abs_error > error_bound) then + has_large_errors = .true. + relative_error = abs_error / max(abs_reference, 1.0e-10) + write(*,*) ' Large error in direction', idir, ' output B(', i, ',', j, '):' + write(*,*) ' Central diff: ', central_diff + write(*,*) ' AD result: ', ad_result + write(*,*) ' Absolute error:', abs_error + write(*,*) ' Error bound:', error_bound + write(*,*) ' Relative error:', relative_error + end if + ! Track max error for reporting (normalized) + relative_error = abs_error / max(abs_reference, 1.0e-10) + max_error = max(max_error, relative_error) + end do + end do + end do + + write(*,*) 'Maximum relative error:', max_error + write(*,*) 'Tolerance thresholds: rtol=1.0e-5, atol=1.0e-5' + if (has_large_errors) then + write(*,*) 'FAIL: Derivatives are outside tolerance' + else + write(*,*) 'PASS: Derivatives are within tolerance (rtol + atol)' + end if + + end subroutine check_derivatives_numerically + +end program test_dtrsm_vector_forward \ No newline at end of file diff --git a/BLAS/test/test_dtrsm_vector_reverse.f90 b/BLAS/test/test_dtrsm_vector_reverse.f90 new file mode 100644 index 0000000..748d90e --- /dev/null +++ b/BLAS/test/test_dtrsm_vector_reverse.f90 @@ -0,0 +1,257 @@ +! Test program for DTRSM vector reverse mode differentiation +! Generated automatically by run_tapenade_blas.py +! Using REAL*8 precision with nbdirsmax=4 + +program test_dtrsm_vector_reverse + implicit none + include 'DIFFSIZES.inc' + + external :: dtrsm + external :: dtrsm_bv + + ! Test parameters + integer, parameter :: n = 4 ! Matrix/vector size for test + integer, parameter :: max_size = n ! Maximum array dimension + integer, parameter :: lda = max_size, ldb = max_size, ldc = max_size ! Leading dimensions + integer :: i, j, k ! Loop counters + integer :: seed_array(33) ! Random seed + real(4) :: temp_real, temp_imag ! Temporary variables for complex initialization + + character :: side + character :: uplo + character :: transa + character :: diag + integer :: msize + integer :: nsize + real(8) :: alpha + real(8), dimension(max_size,max_size) :: a + integer :: lda_val + real(8), dimension(max_size,max_size) :: b + integer :: ldb_val + + ! Adjoint variables (reverse vector mode) + ! In reverse mode: output adjoints are INPUT (cotangents/seeds) + ! input adjoints are OUTPUT (computed gradients) + real(8), dimension(nbdirsmax) :: alphab + real(8), dimension(nbdirsmax,max_size,max_size) :: ab + real(8), dimension(nbdirsmax,max_size,max_size) :: bb + + ! Storage for original cotangents (for INOUT parameters in VJP verification) + real(8), dimension(nbdirsmax,max_size,max_size) :: bb_orig + + ! Storage for original values (for VJP verification) + real(8) :: alpha_orig + real(8), dimension(max_size,max_size) :: a_orig + real(8), dimension(max_size,max_size) :: b_orig + + ! Variables for VJP verification via finite differences + real(8), parameter :: h = 1.0e-7 + real(8) :: vjp_ad, vjp_fd, relative_error, max_error, abs_error, abs_reference, error_bound + logical :: has_large_errors + real(8), dimension(max_size*max_size) :: temp_products ! For sorted summation + integer :: n_products + + ! Initialize random seed for reproducibility + seed_array = 42 + call random_seed(put=seed_array) + + ! Initialize primal values + side = 'L' + uplo = 'U' + transa = 'N' + diag = 'N' + msize = n + nsize = n + call random_number(alpha) + alpha = alpha * 2.0 - 1.0 + call random_number(a) + a = a * 2.0 - 1.0 + lda_val = lda + call random_number(b) + b = b * 2.0 - 1.0 + ldb_val = ldb + + ! Store original primal values + alpha_orig = alpha + a_orig = a + b_orig = b + + ! Initialize output adjoints (cotangents) with random values for each direction + ! These are the 'seeds' for reverse mode + do k = 1, nbdirsmax + call random_number(bb(k,:,:)) + bb(k,:,:) = bb(k,:,:) * 2.0 - 1.0 + end do + + ! Initialize input adjoints to zero (they will be computed) + ! Note: Inout parameters are skipped - they already have output adjoints initialized + alphab = 0.0 + ab = 0.0 + + ! Save original cotangent seeds for OUTPUT/INOUT parameters (before function call) + bb_orig = bb + + ! Set ISIZE globals required by differentiated routine (dimension 2 of arrays). + ! ISIZE1OF* (vectors): use n to match adjoint array size; ISIZE2OF* (matrices): use max_size. + call set_ISIZE2OFA(max_size) + + ! Call reverse vector mode differentiated function + call dtrsm_bv(side, uplo, transa, diag, msize, nsize, alpha, alphab, a, ab, lda_val, b, bb, ldb_val, nbdirsmax) + + ! Reset ISIZE globals to uninitialized (-1) for completeness + call set_ISIZE2OFA(-1) + + ! VJP Verification using finite differences + call check_vjp_numerically() + + write(*,*) '' + write(*,*) 'Test completed successfully' + +contains + + subroutine check_vjp_numerically() + implicit none + + ! Direction vectors for VJP testing + real(8) :: alpha_dir + real(8), dimension(max_size,max_size) :: a_dir + real(8), dimension(max_size,max_size) :: b_dir + real(8), dimension(max_size,max_size) :: b_plus, b_minus, b_central_diff + + max_error = 0.0d0 + has_large_errors = .false. + + write(*,*) 'Function calls completed successfully' + + write(*,*) 'Checking derivatives against numerical differentiation:' + write(*,*) 'Step size h =', h + + ! Test each differentiation direction separately + do k = 1, nbdirsmax + + ! Initialize random direction vectors for all inputs + call random_number(alpha_dir) + alpha_dir = alpha_dir * 2.0 - 1.0 + call random_number(a_dir) + a_dir = a_dir * 2.0 - 1.0 + call random_number(b_dir) + b_dir = b_dir * 2.0 - 1.0 + + ! Forward perturbation: f(x + h*dir) + alpha = alpha_orig + h * alpha_dir + a = a_orig + h * a_dir + b = b_orig + h * b_dir + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + b_plus = b + + ! Backward perturbation: f(x - h*dir) + alpha = alpha_orig - h * alpha_dir + a = a_orig - h * a_dir + b = b_orig - h * b_dir + call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val) + b_minus = b + + ! Compute central differences and VJP verification + ! VJP check: direction^T @ adjoint should equal finite difference + + ! Compute central differences: (f(x+h*dir) - f(x-h*dir)) / (2h) + b_central_diff = (b_plus - b_minus) / (2.0d0 * h) + + ! VJP verification: + ! cotangent^T @ central_diff should equal direction^T @ computed_adjoint + ! Left side: cotangent^T @ Jacobian @ direction (via finite differences, with sorted summation) + vjp_fd = 0.0d0 + ! Compute and sort products for b (FD) + n_products = 0 + do j = 1, n + do i = 1, n + n_products = n_products + 1 + temp_products(n_products) = bb_orig(k,i,j) * b_central_diff(i,j) + end do + end do + call sort_array(temp_products, n_products) + do i = 1, n_products + vjp_fd = vjp_fd + temp_products(i) + end do + + ! Right side: direction^T @ computed_adjoint (with sorted summation) + ! For INOUT parameters: use cb directly (it contains the computed input adjoint after reverse pass) + ! For pure inputs: use adjoint directly + vjp_ad = 0.0d0 + ! Compute and sort products for a + n_products = 0 + do j = 1, n + do i = 1, n + n_products = n_products + 1 + temp_products(n_products) = a_dir(i,j) * ab(k,i,j) + end do + end do + call sort_array(temp_products, n_products) + do i = 1, n_products + vjp_ad = vjp_ad + temp_products(i) + end do + ! Compute and sort products for b + n_products = 0 + do j = 1, n + do i = 1, n + n_products = n_products + 1 + temp_products(n_products) = b_dir(i,j) * bb(k,i,j) + end do + end do + call sort_array(temp_products, n_products) + do i = 1, n_products + vjp_ad = vjp_ad + temp_products(i) + end do + vjp_ad = vjp_ad + alpha_dir * alphab(k) + + ! Error check: |vjp_fd - vjp_ad| > atol + rtol * |vjp_ad| + abs_error = abs(vjp_fd - vjp_ad) + abs_reference = abs(vjp_ad) + error_bound = 1.0e-5 + 1.0e-5 * abs_reference + if (abs_error > error_bound) then + has_large_errors = .true. + end if + + ! Compute relative error for reporting + if (abs_reference > 1.0e-10) then + relative_error = abs_error / abs_reference + else + relative_error = abs_error + end if + if (relative_error > max_error) max_error = relative_error + end do + + write(*,*) 'Maximum relative error:', max_error + write(*,*) 'Tolerance thresholds: rtol=1.0e-5, atol=1.0e-5' + if (has_large_errors) then + write(*,*) 'FAIL: Derivatives are outside tolerance' + else + write(*,*) 'PASS: Derivatives are within tolerance (rtol + atol)' + end if + + end subroutine check_vjp_numerically + + subroutine sort_array(arr, n) + implicit none + integer, intent(in) :: n + real(8), dimension(n), intent(inout) :: arr + integer :: i, j, min_idx + real(8) :: temp + + ! Simple selection sort + do i = 1, n-1 + min_idx = i + do j = i+1, n + if (abs(arr(j)) < abs(arr(min_idx))) then + min_idx = j + end if + end do + if (min_idx /= i) then + temp = arr(i) + arr(i) = arr(min_idx) + arr(min_idx) = temp + end if + end do + end subroutine sort_array + +end program test_dtrsm_vector_reverse diff --git a/README.md b/README.md index 5e5700c..5b145f9 100644 --- a/README.md +++ b/README.md @@ -76,7 +76,7 @@ You need **pre-generated** sources (from step 1). The build compiles them and li **Dependencies:** - **Fortran compiler** (e.g. gfortran, ifort, ifx) and **C compiler** (e.g. gcc). -- **LAPACK installation** — a built Reference LAPACK (or compatible) providing BLAS (e.g. `librefblas.a` or `libblas.a`). Set **`LAPACKDIR`** (or equivalent) so Meson can find it (see below). +- **LAPACK installation** — a built Reference LAPACK (or compatible) providing BLAS (e.g. `librefblas.a` or `libblas.a`). Set **`LAPACKDIR`** (or equivalent) so Meson can find it (see below). If your build produced `libblas.a` rather than `librefblas.a`, pass `-Dlibblas=blas` instead of `-Dlibblas=refblas`. - **Tapenade adStack** — the repo already contains `TAPENADE/adStack.c` and `TAPENADE/include/`; Meson compiles and links these automatically. No separate Tapenade install is required for the build. **Configure and build from the project root:** @@ -116,9 +116,22 @@ meson install -C builddir --prefix /your/install ```bash cd BLAS export LAPACKDIR=/path/to/your/lapack/build # or wherever librefblas is -make +make # use make -j 8 or make -j 16 for a faster build using multi-threading ``` +**Note on library naming:** if you built Reference LAPACK yourself via its standard CMake build, the resulting BLAS archive is usually named `libblas.a` (and `liblapack.a`), not `librefblas.a`. Since this Makefile links against `-lrefblas` by default (see `BLAS_LIB` in the Makefile), you'll need to either: + +1. Symlink it to the expected name: + + cd $LAPACKDIR + ln -s libblas.a librefblas.a + +2. Or override `BLAS_LIB` with the full linker flags for your archive's actual name: + + make BLAS_LIB="-L$LAPACKDIR -lblas" + +Run `find $LAPACKDIR -name "*.a"` first if you're not sure which archive names you actually have. + This builds per-mode static libraries (`build/libdiffblas_d.a`, `libdiffblas_b.a`, `libdiffblas_dv.a`, `libdiffblas_bv.a`) and test executables in `build/`. Run tests: ```bash