diff --git a/lib/linalg/Makefile.gfortran b/lib/linalg/Makefile.gfortran index dca29aa29f..28d1eee07f 100755 --- a/lib/linalg/Makefile.gfortran +++ b/lib/linalg/Makefile.gfortran @@ -1,23 +1,19 @@ # * # *_________________________________________________________________________* -# * Minimal BLAS/LAPACK Library for ATC - -# To compile and link LAMMPS to the linalg library generated by this Makefile, -# try appending the following definitions to the standard definitions in -# whatever LAMMPS Makefile your are using, e.g. in lib/atc/Makefile.lammps - -# CCFLAGS = -I../../lib/linalg -# LINKFLAGS = -L../../lib/libalg -# USRLIB = -llinalg -lgfortran +# * Minimal BLAS/LAPACK Library for use by other LAMMPS packages SHELL = /bin/sh # ------ FILES ------ -SRC = dasum.f daxpy.f dcopy.f ddot.f dgecon.f dgemm.f dgemv.f dger.f dgetf2.f dgetrf.f dgetri.f dlabad.f dlamch.f dlacn2.f dlange.f dlassq.f dlaswp.f dlatrs.f drscl.f dscal.f dswap.f dtrmm.f dtrmv.f dtrsm.f dtrsv.f dtrti2.f dtrtri.f idamax.f ieeeck.f ilaenv.f iparmq.f lsame.f xerbla.f +SRC = dasum.f daxpy.f dcopy.f ddot.f dgecon.f dgemm.f dgemv.f dger.f \ + dgetf2.f dgetrf.f dgetri.f disnan.f dlabad.f dlaisnan.f dlamch.f\ + dlacn2.f dlange.f dlassq.f dlaswp.f dlatrs.f drscl.f dscal.f \ + dswap.f dtrmm.f dtrmv.f dtrsm.f dtrsv.f dtrti2.f dtrtri.f \ + idamax.f ieeeck.f ilaenv.f iparmq.f lsame.f xerbla.f zdotc.f \ + zdscal.f zhpr.f zpptrf.f zpptri.f zscal.f ztpmv.f ztpsv.f ztptri.f - -FILES = $(SRC) Makefile.* +FILES = $(SRC) Makefile.* README # ------ DEFINITIONS ------ @@ -29,7 +25,7 @@ OBJ = $(SRC:.f=.o) FC = gfortran FFLAGS = -O3 -fPIC -march=native -mpc64 \ -ffast-math -funroll-loops -fstrict-aliasing -Wall -W -Wno-uninitialized -fno-second-underscore -FFLAGS0 = -O0 -march=native -mpc64 \ +FFLAGS0 = -O0 -fPIC -march=native -mpc64 \ -Wall -W -Wno-uninitialized -fno-second-underscore ARCHIVE = ar AR = ar diff --git a/lib/linalg/Makefile.mingw_cross b/lib/linalg/Makefile.mingw_cross index 36b5bd5806..40136dbea6 100755 --- a/lib/linalg/Makefile.mingw_cross +++ b/lib/linalg/Makefile.mingw_cross @@ -1,21 +1,17 @@ # * # *_________________________________________________________________________* -# * Minimal BLAS/LAPACK Library for ATC - -# To compile and link LAMMPS to the linalg library generated by this Makefile, -# try appending the following definitions to the standard definitions in -# whatever LAMMPS Makefile your are using, e.g. in lib/atc/Makefile.lammps - -# CCFLAGS = -I../../lib/linalg -# LINKFLAGS = -L../../lib/libalg -# USRLIB = -llinalg -lgfortran +# * Minimal BLAS/LAPACK Library for use by other LAMMPS packages SHELL = /bin/sh # ------ FILES ------ -SRC = dasum.f daxpy.f dcopy.f ddot.f dgecon.f dgemm.f dgemv.f dger.f dgetf2.f dgetrf.f dgetri.f dlabad.f dlamch.f dlacn2.f dlange.f dlassq.f dlaswp.f dlatrs.f drscl.f dscal.f dswap.f dtrmm.f dtrmv.f dtrsm.f dtrsv.f dtrti2.f dtrtri.f idamax.f ieeeck.f ilaenv.f iparmq.f lsame.f xerbla.f - +SRC = dasum.f daxpy.f dcopy.f ddot.f dgecon.f dgemm.f dgemv.f dger.f \ + dgetf2.f dgetrf.f dgetri.f disnan.f dlabad.f dlaisnan.f dlamch.f\ + dlacn2.f dlange.f dlassq.f dlaswp.f dlatrs.f drscl.f dscal.f \ + dswap.f dtrmm.f dtrmv.f dtrsm.f dtrsv.f dtrti2.f dtrtri.f \ + idamax.f ieeeck.f ilaenv.f iparmq.f lsame.f xerbla.f zdotc.f \ + zdscal.f zhpr.f zpptrf.f zpptri.f zscal.f ztpmv.f ztpsv.f ztptri.f FILES = $(SRC) Makefile diff --git a/lib/linalg/README b/lib/linalg/README index 9bfb465f63..20f3ff094d 100644 --- a/lib/linalg/README +++ b/lib/linalg/README @@ -1,9 +1,13 @@ This directory has BLAS and LAPACK files needed by the USER-ATC and USER-AWPMD packages, and possibly by other packages in the future. -You only need to build and use the files in this directory if you want -to build LAMMPS with a package that needs BLAS/LAPACK routines and you -do not have the standard BLAS and LAPACK libraries on your system. +Note that this is an *incomplete* subset of full BLAS/LAPACK. + +You should only need to build and use the resulting library in this +directory if you want to build LAMMPS with the USER-ATC and/or +USER-AWPMD packages AND you do not have any other suitable BLAS and +LAPACK libraries installed on your system. E.g. ATLAS, GOTO-BLAS, +OpenBLAS, ACML, or MKL. Build the library using one of the provided Makefile.* files or create your own, specific to your compiler and system. For example: diff --git a/lib/linalg/dasum.f b/lib/linalg/dasum.f index 8705ecd9a1..c1bd78ac81 100644 --- a/lib/linalg/dasum.f +++ b/lib/linalg/dasum.f @@ -1,4 +1,61 @@ +*> \brief \b DASUM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) +* +* .. Scalar Arguments .. +* INTEGER INCX,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DASUM takes the sum of the absolute values. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 3/93 to return if incx .le. 0. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. INTEGER INCX,N * .. @@ -6,18 +63,6 @@ DOUBLE PRECISION DX(*) * .. * -* Purpose -* ======= -* -* DASUM takes the sum of the absolute values. -* -* Further Details -* =============== -* -* jack dongarra, linpack, 3/11/78. -* modified 3/93 to return if incx .le. 0. -* modified 12/3/93, array(1) declarations changed to array(*) -* * ===================================================================== * * .. Local Scalars .. @@ -30,33 +75,37 @@ DASUM = 0.0d0 DTEMP = 0.0d0 IF (N.LE.0 .OR. INCX.LE.0) RETURN - IF (INCX.EQ.1) GO TO 20 -* -* code for increment not equal to 1 -* - NINCX = N*INCX - DO 10 I = 1,NINCX,INCX - DTEMP = DTEMP + DABS(DX(I)) - 10 CONTINUE - DASUM = DTEMP - RETURN -* + IF (INCX.EQ.1) THEN * code for increment equal to 1 * * * clean-up loop * - 20 M = MOD(N,6) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DTEMP = DTEMP + DABS(DX(I)) - 30 CONTINUE - IF (N.LT.6) GO TO 60 - 40 MP1 = M + 1 - DO 50 I = MP1,N,6 - DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I+1)) + DABS(DX(I+2)) + - + DABS(DX(I+3)) + DABS(DX(I+4)) + DABS(DX(I+5)) - 50 CONTINUE - 60 DASUM = DTEMP + M = MOD(N,6) + IF (M.NE.0) THEN + DO I = 1,M + DTEMP = DTEMP + DABS(DX(I)) + END DO + IF (N.LT.6) THEN + DASUM = DTEMP + RETURN + END IF + END IF + MP1 = M + 1 + DO I = MP1,N,6 + DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I+1)) + + $ DABS(DX(I+2)) + DABS(DX(I+3)) + + $ DABS(DX(I+4)) + DABS(DX(I+5)) + END DO + ELSE +* +* code for increment not equal to 1 +* + NINCX = N*INCX + DO I = 1,NINCX,INCX + DTEMP = DTEMP + DABS(DX(I)) + END DO + END IF + DASUM = DTEMP RETURN END diff --git a/lib/linalg/daxpy.f b/lib/linalg/daxpy.f index 4f92db2fb9..64a02d68bc 100644 --- a/lib/linalg/daxpy.f +++ b/lib/linalg/daxpy.f @@ -1,4 +1,62 @@ +*> \brief \b DAXPY +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION DA +* INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*),DY(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DAXPY constant times a vector plus a vector. +*> uses unrolled loops for increments equal to one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. DOUBLE PRECISION DA INTEGER INCX,INCY,N @@ -7,18 +65,6 @@ DOUBLE PRECISION DX(*),DY(*) * .. * -* Purpose -* ======= -* -* DAXPY constant times a vector plus a vector. -* uses unrolled loops for increments equal to one. -* -* Further Details -* =============== -* -* jack dongarra, linpack, 3/11/78. -* modified 12/3/93, array(1) declarations changed to array(*) -* * ===================================================================== * * .. Local Scalars .. @@ -29,39 +75,41 @@ * .. IF (N.LE.0) RETURN IF (DA.EQ.0.0d0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO 10 I = 1,N - DY(IY) = DY(IY) + DA*DX(IX) - IX = IX + INCX - IY = IY + INCY - 10 CONTINUE - RETURN + IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN * * code for both increments equal to 1 * * * clean-up loop * - 20 M = MOD(N,4) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DY(I) = DY(I) + DA*DX(I) - 30 CONTINUE - IF (N.LT.4) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,4 - DY(I) = DY(I) + DA*DX(I) - DY(I+1) = DY(I+1) + DA*DX(I+1) - DY(I+2) = DY(I+2) + DA*DX(I+2) - DY(I+3) = DY(I+3) + DA*DX(I+3) - 50 CONTINUE + M = MOD(N,4) + IF (M.NE.0) THEN + DO I = 1,M + DY(I) = DY(I) + DA*DX(I) + END DO + END IF + IF (N.LT.4) RETURN + MP1 = M + 1 + DO I = MP1,N,4 + DY(I) = DY(I) + DA*DX(I) + DY(I+1) = DY(I+1) + DA*DX(I+1) + DY(I+2) = DY(I+2) + DA*DX(I+2) + DY(I+3) = DY(I+3) + DA*DX(I+3) + END DO + ELSE +* +* code for unequal increments or equal increments +* not equal to 1 +* + IX = 1 + IY = 1 + IF (INCX.LT.0) IX = (-N+1)*INCX + 1 + IF (INCY.LT.0) IY = (-N+1)*INCY + 1 + DO I = 1,N + DY(IY) = DY(IY) + DA*DX(IX) + IX = IX + INCX + IY = IY + INCY + END DO + END IF RETURN END diff --git a/lib/linalg/dcopy.f b/lib/linalg/dcopy.f index dcc8a0c175..d9d5ac7aa2 100644 --- a/lib/linalg/dcopy.f +++ b/lib/linalg/dcopy.f @@ -1,4 +1,61 @@ +*> \brief \b DCOPY +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) +* +* .. Scalar Arguments .. +* INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*),DY(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DCOPY copies a vector, x, to a vector, y. +*> uses unrolled loops for increments equal to one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. INTEGER INCX,INCY,N * .. @@ -6,18 +63,6 @@ DOUBLE PRECISION DX(*),DY(*) * .. * -* Purpose -* ======= -* -* DCOPY copies a vector, x, to a vector, y. -* uses unrolled loops for increments equal to one. -* -* Further Details -* =============== -* -* jack dongarra, linpack, 3/11/78. -* modified 12/3/93, array(1) declarations changed to array(*) -* * ===================================================================== * * .. Local Scalars .. @@ -27,42 +72,44 @@ INTRINSIC MOD * .. IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO 10 I = 1,N - DY(IY) = DX(IX) - IX = IX + INCX - IY = IY + INCY - 10 CONTINUE - RETURN + IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN * * code for both increments equal to 1 * * * clean-up loop * - 20 M = MOD(N,7) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DY(I) = DX(I) - 30 CONTINUE - IF (N.LT.7) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,7 - DY(I) = DX(I) - DY(I+1) = DX(I+1) - DY(I+2) = DX(I+2) - DY(I+3) = DX(I+3) - DY(I+4) = DX(I+4) - DY(I+5) = DX(I+5) - DY(I+6) = DX(I+6) - 50 CONTINUE + M = MOD(N,7) + IF (M.NE.0) THEN + DO I = 1,M + DY(I) = DX(I) + END DO + IF (N.LT.7) RETURN + END IF + MP1 = M + 1 + DO I = MP1,N,7 + DY(I) = DX(I) + DY(I+1) = DX(I+1) + DY(I+2) = DX(I+2) + DY(I+3) = DX(I+3) + DY(I+4) = DX(I+4) + DY(I+5) = DX(I+5) + DY(I+6) = DX(I+6) + END DO + ELSE +* +* code for unequal increments or equal increments +* not equal to 1 +* + IX = 1 + IY = 1 + IF (INCX.LT.0) IX = (-N+1)*INCX + 1 + IF (INCY.LT.0) IY = (-N+1)*INCY + 1 + DO I = 1,N + DY(IY) = DX(IX) + IX = IX + INCX + IY = IY + INCY + END DO + END IF RETURN END diff --git a/lib/linalg/ddot.f b/lib/linalg/ddot.f index 3df03aecb1..cc0c1b7a43 100644 --- a/lib/linalg/ddot.f +++ b/lib/linalg/ddot.f @@ -1,4 +1,61 @@ +*> \brief \b DDOT +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) +* +* .. Scalar Arguments .. +* INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*),DY(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DDOT forms the dot product of two vectors. +*> uses unrolled loops for increments equal to one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. INTEGER INCX,INCY,N * .. @@ -6,18 +63,6 @@ DOUBLE PRECISION DX(*),DY(*) * .. * -* Purpose -* ======= -* -* DDOT forms the dot product of two vectors. -* uses unrolled loops for increments equal to one. -* -* Further Details -* =============== -* -* jack dongarra, linpack, 3/11/78. -* modified 12/3/93, array(1) declarations changed to array(*) -* * ===================================================================== * * .. Local Scalars .. @@ -30,39 +75,43 @@ DDOT = 0.0d0 DTEMP = 0.0d0 IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 -* -* code for unequal increments or equal increments -* not equal to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO 10 I = 1,N - DTEMP = DTEMP + DX(IX)*DY(IY) - IX = IX + INCX - IY = IY + INCY - 10 CONTINUE - DDOT = DTEMP - RETURN + IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN * * code for both increments equal to 1 * * * clean-up loop * - 20 M = MOD(N,5) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DTEMP = DTEMP + DX(I)*DY(I) - 30 CONTINUE - IF (N.LT.5) GO TO 60 - 40 MP1 = M + 1 - DO 50 I = MP1,N,5 + M = MOD(N,5) + IF (M.NE.0) THEN + DO I = 1,M + DTEMP = DTEMP + DX(I)*DY(I) + END DO + IF (N.LT.5) THEN + DDOT=DTEMP + RETURN + END IF + END IF + MP1 = M + 1 + DO I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) + - + DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) - 50 CONTINUE - 60 DDOT = DTEMP + $ DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) + END DO + ELSE +* +* code for unequal increments or equal increments +* not equal to 1 +* + IX = 1 + IY = 1 + IF (INCX.LT.0) IX = (-N+1)*INCX + 1 + IF (INCY.LT.0) IY = (-N+1)*INCY + 1 + DO I = 1,N + DTEMP = DTEMP + DX(IX)*DY(IY) + IX = IX + INCX + IY = IY + INCY + END DO + END IF + DDOT = DTEMP RETURN END diff --git a/lib/linalg/dgecon.f b/lib/linalg/dgecon.f index b49a91f6f2..df9d8e1c40 100644 --- a/lib/linalg/dgecon.f +++ b/lib/linalg/dgecon.f @@ -1,12 +1,133 @@ +*> \brief \b DGECON +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGECON + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, +* INFO ) +* +* .. Scalar Arguments .. +* CHARACTER NORM +* INTEGER INFO, LDA, N +* DOUBLE PRECISION ANORM, RCOND +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGECON estimates the reciprocal of the condition number of a general +*> real matrix A, in either the 1-norm or the infinity-norm, using +*> the LU factorization computed by DGETRF. +*> +*> An estimate is obtained for norm(inv(A)), and the reciprocal of the +*> condition number is computed as +*> RCOND = 1 / ( norm(A) * norm(inv(A)) ). +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] NORM +*> \verbatim +*> NORM is CHARACTER*1 +*> Specifies whether the 1-norm condition number or the +*> infinity-norm condition number is required: +*> = '1' or 'O': 1-norm; +*> = 'I': Infinity-norm. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> The factors L and U from the factorization A = P*L*U +*> as computed by DGETRF. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] ANORM +*> \verbatim +*> ANORM is DOUBLE PRECISION +*> If NORM = '1' or 'O', the 1-norm of the original matrix A. +*> If NORM = 'I', the infinity-norm of the original matrix A. +*> \endverbatim +*> +*> \param[out] RCOND +*> \verbatim +*> RCOND is DOUBLE PRECISION +*> The reciprocal of the condition number of the matrix A, +*> computed as RCOND = 1/(norm(A) * norm(inv(A))). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (4*N) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleGEcomputational +* +* ===================================================================== SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, $ INFO ) * -* -- LAPACK routine (version 3.2) -- +* -- LAPACK computational routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 -* -* Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. +* November 2011 * * .. Scalar Arguments .. CHARACTER NORM @@ -18,52 +139,6 @@ DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * -* Purpose -* ======= -* -* DGECON estimates the reciprocal of the condition number of a general -* real matrix A, in either the 1-norm or the infinity-norm, using -* the LU factorization computed by DGETRF. -* -* An estimate is obtained for norm(inv(A)), and the reciprocal of the -* condition number is computed as -* RCOND = 1 / ( norm(A) * norm(inv(A)) ). -* -* Arguments -* ========= -* -* NORM (input) CHARACTER*1 -* Specifies whether the 1-norm condition number or the -* infinity-norm condition number is required: -* = '1' or 'O': 1-norm; -* = 'I': Infinity-norm. -* -* N (input) INTEGER -* The order of the matrix A. N >= 0. -* -* A (input) DOUBLE PRECISION array, dimension (LDA,N) -* The factors L and U from the factorization A = P*L*U -* as computed by DGETRF. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* ANORM (input) DOUBLE PRECISION -* If NORM = '1' or 'O', the 1-norm of the original matrix A. -* If NORM = 'I', the infinity-norm of the original matrix A. -* -* RCOND (output) DOUBLE PRECISION -* The reciprocal of the condition number of the matrix A, -* computed as RCOND = 1/(norm(A) * norm(inv(A))). -* -* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) -* -* IWORK (workspace) INTEGER array, dimension (N) -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* * ===================================================================== * * .. Parameters .. @@ -149,12 +224,12 @@ $ A, LDA, WORK, SU, WORK( 3*N+1 ), INFO ) ELSE * -* Multiply by inv(U'). +* Multiply by inv(U**T). * CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A, $ LDA, WORK, SU, WORK( 3*N+1 ), INFO ) * -* Multiply by inv(L'). +* Multiply by inv(L**T). * CALL DLATRS( 'Lower', 'Transpose', 'Unit', NORMIN, N, A, $ LDA, WORK, SL, WORK( 2*N+1 ), INFO ) diff --git a/lib/linalg/dgemm.f b/lib/linalg/dgemm.f index 4afccb1615..45d001b7ab 100644 --- a/lib/linalg/dgemm.f +++ b/lib/linalg/dgemm.f @@ -1,4 +1,197 @@ +*> \brief \b DGEMM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA,BETA +* INTEGER K,LDA,LDB,LDC,M,N +* CHARACTER TRANSA,TRANSB +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGEMM performs one of the matrix-matrix operations +*> +*> C := alpha*op( A )*op( B ) + beta*C, +*> +*> where op( X ) is one of +*> +*> op( X ) = X or op( X ) = X**T, +*> +*> alpha and beta are scalars, and A, B and C are matrices, with op( A ) +*> an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] TRANSA +*> \verbatim +*> TRANSA is CHARACTER*1 +*> On entry, TRANSA specifies the form of op( A ) to be used in +*> the matrix multiplication as follows: +*> +*> TRANSA = 'N' or 'n', op( A ) = A. +*> +*> TRANSA = 'T' or 't', op( A ) = A**T. +*> +*> TRANSA = 'C' or 'c', op( A ) = A**T. +*> \endverbatim +*> +*> \param[in] TRANSB +*> \verbatim +*> TRANSB is CHARACTER*1 +*> On entry, TRANSB specifies the form of op( B ) to be used in +*> the matrix multiplication as follows: +*> +*> TRANSB = 'N' or 'n', op( B ) = B. +*> +*> TRANSB = 'T' or 't', op( B ) = B**T. +*> +*> TRANSB = 'C' or 'c', op( B ) = B**T. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of the matrix +*> op( A ) and of the matrix C. M must be at least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of the matrix +*> op( B ) and the number of columns of the matrix C. N must be +*> at least zero. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> On entry, K specifies the number of columns of the matrix +*> op( A ) and the number of rows of the matrix op( B ). K must +*> be at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is +*> k when TRANSA = 'N' or 'n', and is m otherwise. +*> Before entry with TRANSA = 'N' or 'n', the leading m by k +*> part of the array A must contain the matrix A, otherwise +*> the leading k by m part of the array A must contain the +*> matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. When TRANSA = 'N' or 'n' then +*> LDA must be at least max( 1, m ), otherwise LDA must be at +*> least max( 1, k ). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is +*> n when TRANSB = 'N' or 'n', and is k otherwise. +*> Before entry with TRANSB = 'N' or 'n', the leading k by n +*> part of the array B must contain the matrix B, otherwise +*> the leading n by k part of the array B must contain the +*> matrix B. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> On entry, LDB specifies the first dimension of B as declared +*> in the calling (sub) program. When TRANSB = 'N' or 'n' then +*> LDB must be at least max( 1, k ), otherwise LDB must be at +*> least max( 1, n ). +*> \endverbatim +*> +*> \param[in] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION. +*> On entry, BETA specifies the scalar beta. When BETA is +*> supplied as zero then C need not be set on input. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is DOUBLE PRECISION array of DIMENSION ( LDC, n ). +*> Before entry, the leading m by n part of the array C must +*> contain the matrix C, except when beta is zero, in which +*> case C need not be set on entry. +*> On exit, the array C is overwritten by the m by n matrix +*> ( alpha*op( A )*op( B ) + beta*C ). +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> On entry, LDC specifies the first dimension of C as declared +*> in the calling (sub) program. LDC must be at least +*> max( 1, m ). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level3 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 3 Blas routine. +*> +*> -- Written on 8-February-1989. +*> Jack Dongarra, Argonne National Laboratory. +*> Iain Duff, AERE Harwell. +*> Jeremy Du Croz, Numerical Algorithms Group Ltd. +*> Sven Hammarling, Numerical Algorithms Group Ltd. +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* +* -- Reference BLAS level3 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. DOUBLE PRECISION ALPHA,BETA INTEGER K,LDA,LDB,LDC,M,N @@ -8,127 +201,6 @@ DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) * .. * -* Purpose -* ======= -* -* DGEMM performs one of the matrix-matrix operations -* -* C := alpha*op( A )*op( B ) + beta*C, -* -* where op( X ) is one of -* -* op( X ) = X or op( X ) = X', -* -* alpha and beta are scalars, and A, B and C are matrices, with op( A ) -* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. -* -* Arguments -* ========== -* -* TRANSA - CHARACTER*1. -* On entry, TRANSA specifies the form of op( A ) to be used in -* the matrix multiplication as follows: -* -* TRANSA = 'N' or 'n', op( A ) = A. -* -* TRANSA = 'T' or 't', op( A ) = A'. -* -* TRANSA = 'C' or 'c', op( A ) = A'. -* -* Unchanged on exit. -* -* TRANSB - CHARACTER*1. -* On entry, TRANSB specifies the form of op( B ) to be used in -* the matrix multiplication as follows: -* -* TRANSB = 'N' or 'n', op( B ) = B. -* -* TRANSB = 'T' or 't', op( B ) = B'. -* -* TRANSB = 'C' or 'c', op( B ) = B'. -* -* Unchanged on exit. -* -* M - INTEGER. -* On entry, M specifies the number of rows of the matrix -* op( A ) and of the matrix C. M must be at least zero. -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the number of columns of the matrix -* op( B ) and the number of columns of the matrix C. N must be -* at least zero. -* Unchanged on exit. -* -* K - INTEGER. -* On entry, K specifies the number of columns of the matrix -* op( A ) and the number of rows of the matrix op( B ). K must -* be at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is -* k when TRANSA = 'N' or 'n', and is m otherwise. -* Before entry with TRANSA = 'N' or 'n', the leading m by k -* part of the array A must contain the matrix A, otherwise -* the leading k by m part of the array A must contain the -* matrix A. -* Unchanged on exit. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. When TRANSA = 'N' or 'n' then -* LDA must be at least max( 1, m ), otherwise LDA must be at -* least max( 1, k ). -* Unchanged on exit. -* -* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is -* n when TRANSB = 'N' or 'n', and is k otherwise. -* Before entry with TRANSB = 'N' or 'n', the leading k by n -* part of the array B must contain the matrix B, otherwise -* the leading n by k part of the array B must contain the -* matrix B. -* Unchanged on exit. -* -* LDB - INTEGER. -* On entry, LDB specifies the first dimension of B as declared -* in the calling (sub) program. When TRANSB = 'N' or 'n' then -* LDB must be at least max( 1, k ), otherwise LDB must be at -* least max( 1, n ). -* Unchanged on exit. -* -* BETA - DOUBLE PRECISION. -* On entry, BETA specifies the scalar beta. When BETA is -* supplied as zero then C need not be set on input. -* Unchanged on exit. -* -* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). -* Before entry, the leading m by n part of the array C must -* contain the matrix C, except when beta is zero, in which -* case C need not be set on entry. -* On exit, the array C is overwritten by the m by n matrix -* ( alpha*op( A )*op( B ) + beta*C ). -* -* LDC - INTEGER. -* On entry, LDC specifies the first dimension of C as declared -* in the calling (sub) program. LDC must be at least -* max( 1, m ). -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 3 Blas routine. -* -* -- Written on 8-February-1989. -* Jack Dongarra, Argonne National Laboratory. -* Iain Duff, AERE Harwell. -* Jeremy Du Croz, Numerical Algorithms Group Ltd. -* Sven Hammarling, Numerical Algorithms Group Ltd. -* * ===================================================================== * * .. External Functions .. @@ -249,7 +321,7 @@ 90 CONTINUE ELSE * -* Form C := alpha*A'*B + beta*C +* Form C := alpha*A**T*B + beta*C * DO 120 J = 1,N DO 110 I = 1,M @@ -268,7 +340,7 @@ ELSE IF (NOTA) THEN * -* Form C := alpha*A*B' + beta*C +* Form C := alpha*A*B**T + beta*C * DO 170 J = 1,N IF (BETA.EQ.ZERO) THEN @@ -291,7 +363,7 @@ 170 CONTINUE ELSE * -* Form C := alpha*A'*B' + beta*C +* Form C := alpha*A**T*B**T + beta*C * DO 200 J = 1,N DO 190 I = 1,M diff --git a/lib/linalg/dgemv.f b/lib/linalg/dgemv.f index 59e11b8826..675257fac7 100644 --- a/lib/linalg/dgemv.f +++ b/lib/linalg/dgemv.f @@ -1,4 +1,166 @@ +*> \brief \b DGEMV +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA,BETA +* INTEGER INCX,INCY,LDA,M,N +* CHARACTER TRANS +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),X(*),Y(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGEMV performs one of the matrix-vector operations +*> +*> y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, +*> +*> where alpha and beta are scalars, x and y are vectors and A is an +*> m by n matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> On entry, TRANS specifies the operation to be performed as +*> follows: +*> +*> TRANS = 'N' or 'n' y := alpha*A*x + beta*y. +*> +*> TRANS = 'T' or 't' y := alpha*A**T*x + beta*y. +*> +*> TRANS = 'C' or 'c' y := alpha*A**T*x + beta*y. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of the matrix A. +*> M must be at least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of the matrix A. +*> N must be at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, n ). +*> Before entry, the leading m by n part of the array A must +*> contain the matrix of coefficients. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. LDA must be at least +*> max( 1, m ). +*> \endverbatim +*> +*> \param[in] X +*> \verbatim +*> X is DOUBLE PRECISION array of DIMENSION at least +*> ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' +*> and at least +*> ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. +*> Before entry, the incremented array X must contain the +*> vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> On entry, INCX specifies the increment for the elements of +*> X. INCX must not be zero. +*> \endverbatim +*> +*> \param[in] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION. +*> On entry, BETA specifies the scalar beta. When BETA is +*> supplied as zero then Y need not be set on input. +*> \endverbatim +*> +*> \param[in,out] Y +*> \verbatim +*> Y is DOUBLE PRECISION array of DIMENSION at least +*> ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' +*> and at least +*> ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. +*> Before entry with BETA non-zero, the incremented array Y +*> must contain the vector y. On exit, Y is overwritten by the +*> updated vector y. +*> \endverbatim +*> +*> \param[in] INCY +*> \verbatim +*> INCY is INTEGER +*> On entry, INCY specifies the increment for the elements of +*> Y. INCY must not be zero. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level2 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 2 Blas routine. +*> The vector and matrix arguments are not referenced when N = 0, or M = 0 +*> +*> -- Written on 22-October-1986. +*> Jack Dongarra, Argonne National Lab. +*> Jeremy Du Croz, Nag Central Office. +*> Sven Hammarling, Nag Central Office. +*> Richard Hanson, Sandia National Labs. +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) +* +* -- Reference BLAS level2 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. DOUBLE PRECISION ALPHA,BETA INTEGER INCX,INCY,LDA,M,N @@ -8,98 +170,6 @@ DOUBLE PRECISION A(LDA,*),X(*),Y(*) * .. * -* Purpose -* ======= -* -* DGEMV performs one of the matrix-vector operations -* -* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, -* -* where alpha and beta are scalars, x and y are vectors and A is an -* m by n matrix. -* -* Arguments -* ========== -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the operation to be performed as -* follows: -* -* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. -* -* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. -* -* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. -* -* Unchanged on exit. -* -* M - INTEGER. -* On entry, M specifies the number of rows of the matrix A. -* M must be at least zero. -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the number of columns of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). -* Before entry, the leading m by n part of the array A must -* contain the matrix of coefficients. -* Unchanged on exit. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. LDA must be at least -* max( 1, m ). -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of DIMENSION at least -* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' -* and at least -* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. -* Before entry, the incremented array X must contain the -* vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* BETA - DOUBLE PRECISION. -* On entry, BETA specifies the scalar beta. When BETA is -* supplied as zero then Y need not be set on input. -* Unchanged on exit. -* -* Y - DOUBLE PRECISION array of DIMENSION at least -* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' -* and at least -* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. -* Before entry with BETA non-zero, the incremented array Y -* must contain the vector y. On exit, Y is overwritten by the -* updated vector y. -* -* INCY - INTEGER. -* On entry, INCY specifies the increment for the elements of -* Y. INCY must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* * ===================================================================== * * .. Parameters .. @@ -231,7 +301,7 @@ END IF ELSE * -* Form y := alpha*A'*x + y. +* Form y := alpha*A**T*x + y. * JY = KY IF (INCX.EQ.1) THEN diff --git a/lib/linalg/dger.f b/lib/linalg/dger.f index 0a9dcf2afd..a042483703 100644 --- a/lib/linalg/dger.f +++ b/lib/linalg/dger.f @@ -1,4 +1,140 @@ +*> \brief \b DGER +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA +* INTEGER INCX,INCY,LDA,M,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),X(*),Y(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGER performs the rank 1 operation +*> +*> A := alpha*x*y**T + A, +*> +*> where alpha is a scalar, x is an m element vector, y is an n element +*> vector and A is an m by n matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of the matrix A. +*> M must be at least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of the matrix A. +*> N must be at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. +*> \endverbatim +*> +*> \param[in] X +*> \verbatim +*> X is DOUBLE PRECISION array of dimension at least +*> ( 1 + ( m - 1 )*abs( INCX ) ). +*> Before entry, the incremented array X must contain the m +*> element vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> On entry, INCX specifies the increment for the elements of +*> X. INCX must not be zero. +*> \endverbatim +*> +*> \param[in] Y +*> \verbatim +*> Y is DOUBLE PRECISION array of dimension at least +*> ( 1 + ( n - 1 )*abs( INCY ) ). +*> Before entry, the incremented array Y must contain the n +*> element vector y. +*> \endverbatim +*> +*> \param[in] INCY +*> \verbatim +*> INCY is INTEGER +*> On entry, INCY specifies the increment for the elements of +*> Y. INCY must not be zero. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, n ). +*> Before entry, the leading m by n part of the array A must +*> contain the matrix of coefficients. On exit, A is +*> overwritten by the updated matrix. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. LDA must be at least +*> max( 1, m ). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level2 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 2 Blas routine. +*> +*> -- Written on 22-October-1986. +*> Jack Dongarra, Argonne National Lab. +*> Jeremy Du Croz, Nag Central Office. +*> Sven Hammarling, Nag Central Office. +*> Richard Hanson, Sandia National Labs. +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) +* +* -- Reference BLAS level2 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX,INCY,LDA,M,N @@ -7,77 +143,6 @@ DOUBLE PRECISION A(LDA,*),X(*),Y(*) * .. * -* Purpose -* ======= -* -* DGER performs the rank 1 operation -* -* A := alpha*x*y' + A, -* -* where alpha is a scalar, x is an m element vector, y is an n element -* vector and A is an m by n matrix. -* -* Arguments -* ========== -* -* M - INTEGER. -* On entry, M specifies the number of rows of the matrix A. -* M must be at least zero. -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the number of columns of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of dimension at least -* ( 1 + ( m - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the m -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Y - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCY ) ). -* Before entry, the incremented array Y must contain the n -* element vector y. -* Unchanged on exit. -* -* INCY - INTEGER. -* On entry, INCY specifies the increment for the elements of -* Y. INCY must not be zero. -* Unchanged on exit. -* -* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). -* Before entry, the leading m by n part of the array A must -* contain the matrix of coefficients. On exit, A is -* overwritten by the updated matrix. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. LDA must be at least -* max( 1, m ). -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/dgetf2.f b/lib/linalg/dgetf2.f index a4277c1651..649d0671de 100644 --- a/lib/linalg/dgetf2.f +++ b/lib/linalg/dgetf2.f @@ -1,9 +1,117 @@ +*> \brief \b DGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row interchanges (unblocked algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGETF2 + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGETF2 computes an LU factorization of a general m-by-n matrix A +*> using partial pivoting with row interchanges. +*> +*> The factorization has the form +*> A = P * L * U +*> where P is a permutation matrix, L is lower triangular with unit +*> diagonal elements (lower trapezoidal if m > n), and U is upper +*> triangular (upper trapezoidal if m < n). +*> +*> This is the right-looking Level 2 BLAS version of the algorithm. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the m by n matrix to be factored. +*> On exit, the factors L and U from the factorization +*> A = P*L*U; the unit diagonal elements of L are not stored. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (min(M,N)) +*> The pivot indices; for 1 <= i <= min(M,N), row i of the +*> matrix was interchanged with row IPIV(i). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-th argument had an illegal value +*> > 0: if INFO = k, U(k,k) is exactly zero. The factorization +*> has been completed, but the factor U is exactly +*> singular, and division by zero will occur if it is used +*> to solve a system of equations. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleGEcomputational +* +* ===================================================================== SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) * -* -- LAPACK routine (version 3.2) -- +* -- LAPACK computational routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N @@ -13,49 +121,6 @@ DOUBLE PRECISION A( LDA, * ) * .. * -* Purpose -* ======= -* -* DGETF2 computes an LU factorization of a general m-by-n matrix A -* using partial pivoting with row interchanges. -* -* The factorization has the form -* A = P * L * U -* where P is a permutation matrix, L is lower triangular with unit -* diagonal elements (lower trapezoidal if m > n), and U is upper -* triangular (upper trapezoidal if m < n). -* -* This is the right-looking Level 2 BLAS version of the algorithm. -* -* Arguments -* ========= -* -* M (input) INTEGER -* The number of rows of the matrix A. M >= 0. -* -* N (input) INTEGER -* The number of columns of the matrix A. N >= 0. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) -* On entry, the m by n matrix to be factored. -* On exit, the factors L and U from the factorization -* A = P*L*U; the unit diagonal elements of L are not stored. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,M). -* -* IPIV (output) INTEGER array, dimension (min(M,N)) -* The pivot indices; for 1 <= i <= min(M,N), row i of the -* matrix was interchanged with row IPIV(i). -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -k, the k-th argument had an illegal value -* > 0: if INFO = k, U(k,k) is exactly zero. The factorization -* has been completed, but the factor U is exactly -* singular, and division by zero will occur if it is used -* to solve a system of equations. -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/dgetrf.f b/lib/linalg/dgetrf.f index 3f7aebee98..45bb97f30c 100644 --- a/lib/linalg/dgetrf.f +++ b/lib/linalg/dgetrf.f @@ -1,9 +1,117 @@ +*> \brief \b DGETRF +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGETRF + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGETRF computes an LU factorization of a general M-by-N matrix A +*> using partial pivoting with row interchanges. +*> +*> The factorization has the form +*> A = P * L * U +*> where P is a permutation matrix, L is lower triangular with unit +*> diagonal elements (lower trapezoidal if m > n), and U is upper +*> triangular (upper trapezoidal if m < n). +*> +*> This is the right-looking Level 3 BLAS version of the algorithm. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the M-by-N matrix to be factored. +*> On exit, the factors L and U from the factorization +*> A = P*L*U; the unit diagonal elements of L are not stored. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (min(M,N)) +*> The pivot indices; for 1 <= i <= min(M,N), row i of the +*> matrix was interchanged with row IPIV(i). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, U(i,i) is exactly zero. The factorization +*> has been completed, but the factor U is exactly +*> singular, and division by zero will occur if it is used +*> to solve a system of equations. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleGEcomputational +* +* ===================================================================== SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) * -* -- LAPACK routine (version 3.2) -- +* -- LAPACK computational routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N @@ -13,49 +121,6 @@ DOUBLE PRECISION A( LDA, * ) * .. * -* Purpose -* ======= -* -* DGETRF computes an LU factorization of a general M-by-N matrix A -* using partial pivoting with row interchanges. -* -* The factorization has the form -* A = P * L * U -* where P is a permutation matrix, L is lower triangular with unit -* diagonal elements (lower trapezoidal if m > n), and U is upper -* triangular (upper trapezoidal if m < n). -* -* This is the right-looking Level 3 BLAS version of the algorithm. -* -* Arguments -* ========= -* -* M (input) INTEGER -* The number of rows of the matrix A. M >= 0. -* -* N (input) INTEGER -* The number of columns of the matrix A. N >= 0. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) -* On entry, the M-by-N matrix to be factored. -* On exit, the factors L and U from the factorization -* A = P*L*U; the unit diagonal elements of L are not stored. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,M). -* -* IPIV (output) INTEGER array, dimension (min(M,N)) -* The pivot indices; for 1 <= i <= min(M,N), row i of the -* matrix was interchanged with row IPIV(i). -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* > 0: if INFO = i, U(i,i) is exactly zero. The factorization -* has been completed, but the factor U is exactly -* singular, and division by zero will occur if it is used -* to solve a system of equations. -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/dgetri.f b/lib/linalg/dgetri.f index 09367b64e1..ad5324c07e 100644 --- a/lib/linalg/dgetri.f +++ b/lib/linalg/dgetri.f @@ -1,9 +1,123 @@ +*> \brief \b DGETRI +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGETRI + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LDA, LWORK, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGETRI computes the inverse of a matrix using the LU factorization +*> computed by DGETRF. +*> +*> This method inverts U and then computes inv(A) by solving the system +*> inv(A)*L = inv(U) for inv(A). +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the factors L and U from the factorization +*> A = P*L*U as computed by DGETRF. +*> On exit, if INFO = 0, the inverse of the original matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices from DGETRF; for 1<=i<=N, row i of the +*> matrix was interchanged with row IPIV(i). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +*> On exit, if INFO=0, then WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK. LWORK >= max(1,N). +*> For optimal performance LWORK >= N*NB, where NB is +*> the optimal blocksize returned by ILAENV. +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, U(i,i) is exactly zero; the matrix is +*> singular and its inverse could not be computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleGEcomputational +* +* ===================================================================== SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) * -* -- LAPACK routine (version 3.2) -- +* -- LAPACK computational routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, N @@ -13,52 +127,6 @@ DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * -* Purpose -* ======= -* -* DGETRI computes the inverse of a matrix using the LU factorization -* computed by DGETRF. -* -* This method inverts U and then computes inv(A) by solving the system -* inv(A)*L = inv(U) for inv(A). -* -* Arguments -* ========= -* -* N (input) INTEGER -* The order of the matrix A. N >= 0. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) -* On entry, the factors L and U from the factorization -* A = P*L*U as computed by DGETRF. -* On exit, if INFO = 0, the inverse of the original matrix A. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* IPIV (input) INTEGER array, dimension (N) -* The pivot indices from DGETRF; for 1<=i<=N, row i of the -* matrix was interchanged with row IPIV(i). -* -* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) -* On exit, if INFO=0, then WORK(1) returns the optimal LWORK. -* -* LWORK (input) INTEGER -* The dimension of the array WORK. LWORK >= max(1,N). -* For optimal performance LWORK >= N*NB, where NB is -* the optimal blocksize returned by ILAENV. -* -* If LWORK = -1, then a workspace query is assumed; the routine -* only calculates the optimal size of the WORK array, returns -* this value as the first entry of the WORK array, and no error -* message related to LWORK is issued by XERBLA. -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is -* singular and its inverse could not be computed. -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/dlabad.f b/lib/linalg/dlabad.f index 8eb5ce4fcf..9eda3c91db 100644 --- a/lib/linalg/dlabad.f +++ b/lib/linalg/dlabad.f @@ -1,39 +1,88 @@ +*> \brief \b DLABAD +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLABAD + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLABAD( SMALL, LARGE ) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION LARGE, SMALL +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLABAD takes as input the values computed by DLAMCH for underflow and +*> overflow, and returns the square root of each of these values if the +*> log of LARGE is sufficiently large. This subroutine is intended to +*> identify machines with a large exponent range, such as the Crays, and +*> redefine the underflow and overflow limits to be the square roots of +*> the values computed by DLAMCH. This subroutine is needed because +*> DLAMCH does not compensate for poor arithmetic in the upper half of +*> the exponent range, as is found on a Cray. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in,out] SMALL +*> \verbatim +*> SMALL is DOUBLE PRECISION +*> On entry, the underflow threshold as computed by DLAMCH. +*> On exit, if LOG10(LARGE) is sufficiently large, the square +*> root of SMALL, otherwise unchanged. +*> \endverbatim +*> +*> \param[in,out] LARGE +*> \verbatim +*> LARGE is DOUBLE PRECISION +*> On entry, the overflow threshold as computed by DLAMCH. +*> On exit, if LOG10(LARGE) is sufficiently large, the square +*> root of LARGE, otherwise unchanged. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== SUBROUTINE DLABAD( SMALL, LARGE ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. DOUBLE PRECISION LARGE, SMALL * .. * -* Purpose -* ======= -* -* DLABAD takes as input the values computed by DLAMCH for underflow and -* overflow, and returns the square root of each of these values if the -* log of LARGE is sufficiently large. This subroutine is intended to -* identify machines with a large exponent range, such as the Crays, and -* redefine the underflow and overflow limits to be the square roots of -* the values computed by DLAMCH. This subroutine is needed because -* DLAMCH does not compensate for poor arithmetic in the upper half of -* the exponent range, as is found on a Cray. -* -* Arguments -* ========= -* -* SMALL (input/output) DOUBLE PRECISION -* On entry, the underflow threshold as computed by DLAMCH. -* On exit, if LOG10(LARGE) is sufficiently large, the square -* root of SMALL, otherwise unchanged. -* -* LARGE (input/output) DOUBLE PRECISION -* On entry, the overflow threshold as computed by DLAMCH. -* On exit, if LOG10(LARGE) is sufficiently large, the square -* root of LARGE, otherwise unchanged. -* * ===================================================================== * * .. Intrinsic Functions .. diff --git a/lib/linalg/dlacn2.f b/lib/linalg/dlacn2.f index 26b8f3c75b..9dd3c85ea2 100644 --- a/lib/linalg/dlacn2.f +++ b/lib/linalg/dlacn2.f @@ -1,9 +1,145 @@ +*> \brief \b DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLACN2 + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE ) +* +* .. Scalar Arguments .. +* INTEGER KASE, N +* DOUBLE PRECISION EST +* .. +* .. Array Arguments .. +* INTEGER ISGN( * ), ISAVE( 3 ) +* DOUBLE PRECISION V( * ), X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLACN2 estimates the 1-norm of a square, real matrix A. +*> Reverse communication is used for evaluating matrix-vector products. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix. N >= 1. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension (N) +*> On the final return, V = A*W, where EST = norm(V)/norm(W) +*> (W is not returned). +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is DOUBLE PRECISION array, dimension (N) +*> On an intermediate return, X should be overwritten by +*> A * X, if KASE=1, +*> A**T * X, if KASE=2, +*> and DLACN2 must be re-called with all the other parameters +*> unchanged. +*> \endverbatim +*> +*> \param[out] ISGN +*> \verbatim +*> ISGN is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[in,out] EST +*> \verbatim +*> EST is DOUBLE PRECISION +*> On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be +*> unchanged from the previous call to DLACN2. +*> On exit, EST is an estimate (a lower bound) for norm(A). +*> \endverbatim +*> +*> \param[in,out] KASE +*> \verbatim +*> KASE is INTEGER +*> On the initial call to DLACN2, KASE should be 0. +*> On an intermediate return, KASE will be 1 or 2, indicating +*> whether X should be overwritten by A * X or A**T * X. +*> On the final return from DLACN2, KASE will again be 0. +*> \endverbatim +*> +*> \param[in,out] ISAVE +*> \verbatim +*> ISAVE is INTEGER array, dimension (3) +*> ISAVE is used to save variables between calls to DLACN2 +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Originally named SONEST, dated March 16, 1988. +*> +*> This is a thread safe version of DLACON, which uses the array ISAVE +*> in place of a SAVE statement, as follows: +*> +*> DLACON DLACN2 +*> JUMP ISAVE(1) +*> J ISAVE(2) +*> ITER ISAVE(3) +*> \endverbatim +* +*> \par Contributors: +* ================== +*> +*> Nick Higham, University of Manchester +* +*> \par References: +* ================ +*> +*> N.J. Higham, "FORTRAN codes for estimating the one-norm of +*> a real or complex matrix, with applications to condition estimation", +*> ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. +*> +* ===================================================================== SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. INTEGER KASE, N @@ -14,63 +150,6 @@ DOUBLE PRECISION V( * ), X( * ) * .. * -* Purpose -* ======= -* -* DLACN2 estimates the 1-norm of a square, real matrix A. -* Reverse communication is used for evaluating matrix-vector products. -* -* Arguments -* ========= -* -* N (input) INTEGER -* The order of the matrix. N >= 1. -* -* V (workspace) DOUBLE PRECISION array, dimension (N) -* On the final return, V = A*W, where EST = norm(V)/norm(W) -* (W is not returned). -* -* X (input/output) DOUBLE PRECISION array, dimension (N) -* On an intermediate return, X should be overwritten by -* A * X, if KASE=1, -* A' * X, if KASE=2, -* and DLACN2 must be re-called with all the other parameters -* unchanged. -* -* ISGN (workspace) INTEGER array, dimension (N) -* -* EST (input/output) DOUBLE PRECISION -* On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be -* unchanged from the previous call to DLACN2. -* On exit, EST is an estimate (a lower bound) for norm(A). -* -* KASE (input/output) INTEGER -* On the initial call to DLACN2, KASE should be 0. -* On an intermediate return, KASE will be 1 or 2, indicating -* whether X should be overwritten by A * X or A' * X. -* On the final return from DLACN2, KASE will again be 0. -* -* ISAVE (input/output) INTEGER array, dimension (3) -* ISAVE is used to save variables between calls to DLACN2 -* -* Further Details -* ======= ======= -* -* Contributed by Nick Higham, University of Manchester. -* Originally named SONEST, dated March 16, 1988. -* -* Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of -* a real or complex matrix, with applications to condition estimation", -* ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. -* -* This is a thread safe version of DLACON, which uses the array ISAVE -* in place of a SAVE statement, as follows: -* -* DLACON DLACN2 -* JUMP ISAVE(1) -* J ISAVE(2) -* ITER ISAVE(3) -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/dlange.f b/lib/linalg/dlange.f index 471c2371c7..bec815d1ef 100644 --- a/lib/linalg/dlange.f +++ b/lib/linalg/dlange.f @@ -1,9 +1,123 @@ +*> \brief \b DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general rectangular matrix. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLANGE + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* DOUBLE PRECISION FUNCTION DLANGE( NORM, M, N, A, LDA, WORK ) +* +* .. Scalar Arguments .. +* CHARACTER NORM +* INTEGER LDA, M, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLANGE returns the value of the one norm, or the Frobenius norm, or +*> the infinity norm, or the element of largest absolute value of a +*> real matrix A. +*> \endverbatim +*> +*> \return DLANGE +*> \verbatim +*> +*> DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm' +*> ( +*> ( norm1(A), NORM = '1', 'O' or 'o' +*> ( +*> ( normI(A), NORM = 'I' or 'i' +*> ( +*> ( normF(A), NORM = 'F', 'f', 'E' or 'e' +*> +*> where norm1 denotes the one norm of a matrix (maximum column sum), +*> normI denotes the infinity norm of a matrix (maximum row sum) and +*> normF denotes the Frobenius norm of a matrix (square root of sum of +*> squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] NORM +*> \verbatim +*> NORM is CHARACTER*1 +*> Specifies the value to be returned in DLANGE as described +*> above. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. When M = 0, +*> DLANGE is set to zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix A. N >= 0. When N = 0, +*> DLANGE is set to zero. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> The m by n matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(M,1). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), +*> where LWORK >= M when NORM = 'I'; otherwise, WORK is not +*> referenced. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleGEauxiliary +* +* ===================================================================== DOUBLE PRECISION FUNCTION DLANGE( NORM, M, N, A, LDA, WORK ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. CHARACTER NORM @@ -13,56 +127,6 @@ DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * -* Purpose -* ======= -* -* DLANGE returns the value of the one norm, or the Frobenius norm, or -* the infinity norm, or the element of largest absolute value of a -* real matrix A. -* -* Description -* =========== -* -* DLANGE returns the value -* -* DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm' -* ( -* ( norm1(A), NORM = '1', 'O' or 'o' -* ( -* ( normI(A), NORM = 'I' or 'i' -* ( -* ( normF(A), NORM = 'F', 'f', 'E' or 'e' -* -* where norm1 denotes the one norm of a matrix (maximum column sum), -* normI denotes the infinity norm of a matrix (maximum row sum) and -* normF denotes the Frobenius norm of a matrix (square root of sum of -* squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. -* -* Arguments -* ========= -* -* NORM (input) CHARACTER*1 -* Specifies the value to be returned in DLANGE as described -* above. -* -* M (input) INTEGER -* The number of rows of the matrix A. M >= 0. When M = 0, -* DLANGE is set to zero. -* -* N (input) INTEGER -* The number of columns of the matrix A. N >= 0. When N = 0, -* DLANGE is set to zero. -* -* A (input) DOUBLE PRECISION array, dimension (LDA,N) -* The m by n matrix A. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(M,1). -* -* WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)), -* where LWORK >= M when NORM = 'I'; otherwise, WORK is not -* referenced. -* * ===================================================================== * * .. Parameters .. @@ -71,17 +135,17 @@ * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SCALE, SUM, VALUE, TEMP * .. * .. External Subroutines .. EXTERNAL DLASSQ * .. * .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME + LOGICAL LSAME, DISNAN + EXTERNAL LSAME, DISNAN * .. * .. Intrinsic Functions .. - INTRINSIC ABS, MAX, MIN, SQRT + INTRINSIC ABS, MIN, SQRT * .. * .. Executable Statements .. * @@ -94,7 +158,8 @@ VALUE = ZERO DO 20 J = 1, N DO 10 I = 1, M - VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + TEMP = ABS( A( I, J ) ) + IF( VALUE.LT.TEMP .OR. DISNAN( TEMP ) ) VALUE = TEMP 10 CONTINUE 20 CONTINUE ELSE IF( ( LSAME( NORM, 'O' ) ) .OR. ( NORM.EQ.'1' ) ) THEN @@ -107,7 +172,7 @@ DO 30 I = 1, M SUM = SUM + ABS( A( I, J ) ) 30 CONTINUE - VALUE = MAX( VALUE, SUM ) + IF( VALUE.LT.SUM .OR. DISNAN( SUM ) ) VALUE = SUM 40 CONTINUE ELSE IF( LSAME( NORM, 'I' ) ) THEN * @@ -123,7 +188,8 @@ 70 CONTINUE VALUE = ZERO DO 80 I = 1, M - VALUE = MAX( VALUE, WORK( I ) ) + TEMP = WORK( I ) + IF( VALUE.LT.TEMP .OR. DISNAN( TEMP ) ) VALUE = TEMP 80 CONTINUE ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * diff --git a/lib/linalg/dlassq.f b/lib/linalg/dlassq.f index 51cce918df..c7c4087e80 100644 --- a/lib/linalg/dlassq.f +++ b/lib/linalg/dlassq.f @@ -1,9 +1,112 @@ +*> \brief \b DLASSQ updates a sum of squares represented in scaled form. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLASSQ + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* DOUBLE PRECISION SCALE, SUMSQ +* .. +* .. Array Arguments .. +* DOUBLE PRECISION X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLASSQ returns the values scl and smsq such that +*> +*> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +*> +*> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +*> assumed to be non-negative and scl returns the value +*> +*> scl = max( scale, abs( x( i ) ) ). +*> +*> scale and sumsq must be supplied in SCALE and SUMSQ and +*> scl and smsq are overwritten on SCALE and SUMSQ respectively. +*> +*> The routine makes only one pass through the vector x. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of elements to be used from the vector X. +*> \endverbatim +*> +*> \param[in] X +*> \verbatim +*> X is DOUBLE PRECISION array, dimension (N) +*> The vector for which a scaled sum of squares is computed. +*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between successive values of the vector X. +*> INCX > 0. +*> \endverbatim +*> +*> \param[in,out] SCALE +*> \verbatim +*> SCALE is DOUBLE PRECISION +*> On entry, the value scale in the equation above. +*> On exit, SCALE is overwritten with scl , the scaling factor +*> for the sum of squares. +*> \endverbatim +*> +*> \param[in,out] SUMSQ +*> \verbatim +*> SUMSQ is DOUBLE PRECISION +*> On entry, the value sumsq in the equation above. +*> On exit, SUMSQ is overwritten with smsq , the basic sum of +*> squares from which scl has been factored out. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. INTEGER INCX, N @@ -13,47 +116,6 @@ DOUBLE PRECISION X( * ) * .. * -* Purpose -* ======= -* -* DLASSQ returns the values scl and smsq such that -* -* ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -* -* where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is -* assumed to be non-negative and scl returns the value -* -* scl = max( scale, abs( x( i ) ) ). -* -* scale and sumsq must be supplied in SCALE and SUMSQ and -* scl and smsq are overwritten on SCALE and SUMSQ respectively. -* -* The routine makes only one pass through the vector x. -* -* Arguments -* ========= -* -* N (input) INTEGER -* The number of elements to be used from the vector X. -* -* X (input) DOUBLE PRECISION array, dimension (N) -* The vector for which a scaled sum of squares is computed. -* x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. -* -* INCX (input) INTEGER -* The increment between successive values of the vector X. -* INCX > 0. -* -* SCALE (input/output) DOUBLE PRECISION -* On entry, the value scale in the equation above. -* On exit, SCALE is overwritten with scl , the scaling factor -* for the sum of squares. -* -* SUMSQ (input/output) DOUBLE PRECISION -* On entry, the value sumsq in the equation above. -* On exit, SUMSQ is overwritten with smsq , the basic sum of -* squares from which scl has been factored out. -* * ===================================================================== * * .. Parameters .. @@ -64,6 +126,10 @@ INTEGER IX DOUBLE PRECISION ABSXI * .. +* .. External Functions .. + LOGICAL DISNAN + EXTERNAL DISNAN +* .. * .. Intrinsic Functions .. INTRINSIC ABS * .. @@ -71,8 +137,8 @@ * IF( N.GT.0 ) THEN DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX - IF( X( IX ).NE.ZERO ) THEN - ABSXI = ABS( X( IX ) ) + ABSXI = ABS( X( IX ) ) + IF( ABSXI.GT.ZERO.OR.DISNAN( ABSXI ) ) THEN IF( SCALE.LT.ABSXI ) THEN SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2 SCALE = ABSXI diff --git a/lib/linalg/dlaswp.f b/lib/linalg/dlaswp.f index 5239ba1f5f..937e12b2f0 100644 --- a/lib/linalg/dlaswp.f +++ b/lib/linalg/dlaswp.f @@ -1,9 +1,123 @@ +*> \brief \b DLASWP performs a series of row interchanges on a general rectangular matrix. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLASWP + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX ) +* +* .. Scalar Arguments .. +* INTEGER INCX, K1, K2, LDA, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLASWP performs a series of row interchanges on the matrix A. +*> One row interchange is initiated for each of rows K1 through K2 of A. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix A. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the matrix of column dimension N to which the row +*> interchanges will be applied. +*> On exit, the permuted matrix. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. +*> \endverbatim +*> +*> \param[in] K1 +*> \verbatim +*> K1 is INTEGER +*> The first element of IPIV for which a row interchange will +*> be done. +*> \endverbatim +*> +*> \param[in] K2 +*> \verbatim +*> K2 is INTEGER +*> The last element of IPIV for which a row interchange will +*> be done. +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (K2*abs(INCX)) +*> The vector of pivot indices. Only the elements in positions +*> K1 through K2 of IPIV are accessed. +*> IPIV(K) = L implies rows K and L are to be interchanged. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between successive values of IPIV. If IPIV +*> is negative, the pivots are applied in reverse order. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Modified by +*> R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. INTEGER INCX, K1, K2, LDA, N @@ -13,49 +127,6 @@ DOUBLE PRECISION A( LDA, * ) * .. * -* Purpose -* ======= -* -* DLASWP performs a series of row interchanges on the matrix A. -* One row interchange is initiated for each of rows K1 through K2 of A. -* -* Arguments -* ========= -* -* N (input) INTEGER -* The number of columns of the matrix A. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) -* On entry, the matrix of column dimension N to which the row -* interchanges will be applied. -* On exit, the permuted matrix. -* -* LDA (input) INTEGER -* The leading dimension of the array A. -* -* K1 (input) INTEGER -* The first element of IPIV for which a row interchange will -* be done. -* -* K2 (input) INTEGER -* The last element of IPIV for which a row interchange will -* be done. -* -* IPIV (input) INTEGER array, dimension (K2*abs(INCX)) -* The vector of pivot indices. Only the elements in positions -* K1 through K2 of IPIV are accessed. -* IPIV(K) = L implies rows K and L are to be interchanged. -* -* INCX (input) INTEGER -* The increment between successive values of IPIV. If IPIV -* is negative, the pivots are applied in reverse order. -* -* Further Details -* =============== -* -* Modified by -* R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA -* * ===================================================================== * * .. Local Scalars .. diff --git a/lib/linalg/dlatrs.f b/lib/linalg/dlatrs.f index de70c1f565..b34795eb15 100644 --- a/lib/linalg/dlatrs.f +++ b/lib/linalg/dlatrs.f @@ -1,10 +1,247 @@ +*> \brief \b DLATRS solves a triangular system of equations with the scale factor set to prevent overflow. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLATRS + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, +* CNORM, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER DIAG, NORMIN, TRANS, UPLO +* INTEGER INFO, LDA, N +* DOUBLE PRECISION SCALE +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLATRS solves one of the triangular systems +*> +*> A *x = s*b or A**T *x = s*b +*> +*> with scaling to prevent overflow. Here A is an upper or lower +*> triangular matrix, A**T denotes the transpose of A, x and b are +*> n-element vectors, and s is a scaling factor, usually less than +*> or equal to 1, chosen so that the components of x will be less than +*> the overflow threshold. If the unscaled problem will not cause +*> overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A +*> is singular (A(j,j) = 0 for some j), then s is set to 0 and a +*> non-trivial solution to A*x = 0 is returned. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the matrix A is upper or lower triangular. +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> Specifies the operation applied to A. +*> = 'N': Solve A * x = s*b (No transpose) +*> = 'T': Solve A**T* x = s*b (Transpose) +*> = 'C': Solve A**T* x = s*b (Conjugate transpose = Transpose) +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> Specifies whether or not the matrix A is unit triangular. +*> = 'N': Non-unit triangular +*> = 'U': Unit triangular +*> \endverbatim +*> +*> \param[in] NORMIN +*> \verbatim +*> NORMIN is CHARACTER*1 +*> Specifies whether CNORM has been set or not. +*> = 'Y': CNORM contains the column norms on entry +*> = 'N': CNORM is not set on entry. On exit, the norms will +*> be computed and stored in CNORM. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> The triangular matrix A. If UPLO = 'U', the leading n by n +*> upper triangular part of the array A contains the upper +*> triangular matrix, and the strictly lower triangular part of +*> A is not referenced. If UPLO = 'L', the leading n by n lower +*> triangular part of the array A contains the lower triangular +*> matrix, and the strictly upper triangular part of A is not +*> referenced. If DIAG = 'U', the diagonal elements of A are +*> also not referenced and are assumed to be 1. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max (1,N). +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is DOUBLE PRECISION array, dimension (N) +*> On entry, the right hand side b of the triangular system. +*> On exit, X is overwritten by the solution vector x. +*> \endverbatim +*> +*> \param[out] SCALE +*> \verbatim +*> SCALE is DOUBLE PRECISION +*> The scaling factor s for the triangular system +*> A * x = s*b or A**T* x = s*b. +*> If SCALE = 0, the matrix A is singular or badly scaled, and +*> the vector x is an exact or approximate solution to A*x = 0. +*> \endverbatim +*> +*> \param[in,out] CNORM +*> \verbatim +*> CNORM is DOUBLE PRECISION array, dimension (N) +*> +*> If NORMIN = 'Y', CNORM is an input argument and CNORM(j) +*> contains the norm of the off-diagonal part of the j-th column +*> of A. If TRANS = 'N', CNORM(j) must be greater than or equal +*> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) +*> must be greater than or equal to the 1-norm. +*> +*> If NORMIN = 'N', CNORM is an output argument and CNORM(j) +*> returns the 1-norm of the offdiagonal part of the j-th column +*> of A. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-th argument had an illegal value +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> A rough bound on x is computed; if that is less than overflow, DTRSV +*> is called, otherwise, specific code is used which checks for possible +*> overflow or divide-by-zero at every operation. +*> +*> A columnwise scheme is used for solving A*x = b. The basic algorithm +*> if A is lower triangular is +*> +*> x[1:n] := b[1:n] +*> for j = 1, ..., n +*> x(j) := x(j) / A(j,j) +*> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] +*> end +*> +*> Define bounds on the components of x after j iterations of the loop: +*> M(j) = bound on x[1:j] +*> G(j) = bound on x[j+1:n] +*> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. +*> +*> Then for iteration j+1 we have +*> M(j+1) <= G(j) / | A(j+1,j+1) | +*> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | +*> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) +*> +*> where CNORM(j+1) is greater than or equal to the infinity-norm of +*> column j+1 of A, not counting the diagonal. Hence +*> +*> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) +*> 1<=i<=j +*> and +*> +*> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) +*> 1<=i< j +*> +*> Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the +*> reciprocal of the largest M(j), j=1,..,n, is larger than +*> max(underflow, 1/overflow). +*> +*> The bound on x(j) is also used to determine when a step in the +*> columnwise method can be performed without fear of overflow. If +*> the computed bound is greater than a large constant, x is scaled to +*> prevent overflow, but if the bound overflows, x is set to 0, x(j) to +*> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. +*> +*> Similarly, a row-wise scheme is used to solve A**T*x = b. The basic +*> algorithm for A upper triangular is +*> +*> for j = 1, ..., n +*> x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j) +*> end +*> +*> We simultaneously compute two bounds +*> G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j +*> M(j) = bound on x(i), 1<=i<=j +*> +*> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we +*> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. +*> Then the bound on x(j) is +*> +*> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | +*> +*> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) +*> 1<=i<=j +*> +*> and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater +*> than max(underflow, 1/overflow). +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, $ CNORM, INFO ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. CHARACTER DIAG, NORMIN, TRANS, UPLO @@ -15,158 +252,6 @@ DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * ) * .. * -* Purpose -* ======= -* -* DLATRS solves one of the triangular systems -* -* A *x = s*b or A'*x = s*b -* -* with scaling to prevent overflow. Here A is an upper or lower -* triangular matrix, A' denotes the transpose of A, x and b are -* n-element vectors, and s is a scaling factor, usually less than -* or equal to 1, chosen so that the components of x will be less than -* the overflow threshold. If the unscaled problem will not cause -* overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A -* is singular (A(j,j) = 0 for some j), then s is set to 0 and a -* non-trivial solution to A*x = 0 is returned. -* -* Arguments -* ========= -* -* UPLO (input) CHARACTER*1 -* Specifies whether the matrix A is upper or lower triangular. -* = 'U': Upper triangular -* = 'L': Lower triangular -* -* TRANS (input) CHARACTER*1 -* Specifies the operation applied to A. -* = 'N': Solve A * x = s*b (No transpose) -* = 'T': Solve A'* x = s*b (Transpose) -* = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose) -* -* DIAG (input) CHARACTER*1 -* Specifies whether or not the matrix A is unit triangular. -* = 'N': Non-unit triangular -* = 'U': Unit triangular -* -* NORMIN (input) CHARACTER*1 -* Specifies whether CNORM has been set or not. -* = 'Y': CNORM contains the column norms on entry -* = 'N': CNORM is not set on entry. On exit, the norms will -* be computed and stored in CNORM. -* -* N (input) INTEGER -* The order of the matrix A. N >= 0. -* -* A (input) DOUBLE PRECISION array, dimension (LDA,N) -* The triangular matrix A. If UPLO = 'U', the leading n by n -* upper triangular part of the array A contains the upper -* triangular matrix, and the strictly lower triangular part of -* A is not referenced. If UPLO = 'L', the leading n by n lower -* triangular part of the array A contains the lower triangular -* matrix, and the strictly upper triangular part of A is not -* referenced. If DIAG = 'U', the diagonal elements of A are -* also not referenced and are assumed to be 1. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max (1,N). -* -* X (input/output) DOUBLE PRECISION array, dimension (N) -* On entry, the right hand side b of the triangular system. -* On exit, X is overwritten by the solution vector x. -* -* SCALE (output) DOUBLE PRECISION -* The scaling factor s for the triangular system -* A * x = s*b or A'* x = s*b. -* If SCALE = 0, the matrix A is singular or badly scaled, and -* the vector x is an exact or approximate solution to A*x = 0. -* -* CNORM (input or output) DOUBLE PRECISION array, dimension (N) -* -* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) -* contains the norm of the off-diagonal part of the j-th column -* of A. If TRANS = 'N', CNORM(j) must be greater than or equal -* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) -* must be greater than or equal to the 1-norm. -* -* If NORMIN = 'N', CNORM is an output argument and CNORM(j) -* returns the 1-norm of the offdiagonal part of the j-th column -* of A. -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -k, the k-th argument had an illegal value -* -* Further Details -* ======= ======= -* -* A rough bound on x is computed; if that is less than overflow, DTRSV -* is called, otherwise, specific code is used which checks for possible -* overflow or divide-by-zero at every operation. -* -* A columnwise scheme is used for solving A*x = b. The basic algorithm -* if A is lower triangular is -* -* x[1:n] := b[1:n] -* for j = 1, ..., n -* x(j) := x(j) / A(j,j) -* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] -* end -* -* Define bounds on the components of x after j iterations of the loop: -* M(j) = bound on x[1:j] -* G(j) = bound on x[j+1:n] -* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. -* -* Then for iteration j+1 we have -* M(j+1) <= G(j) / | A(j+1,j+1) | -* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | -* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) -* -* where CNORM(j+1) is greater than or equal to the infinity-norm of -* column j+1 of A, not counting the diagonal. Hence -* -* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) -* 1<=i<=j -* and -* -* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) -* 1<=i< j -* -* Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the -* reciprocal of the largest M(j), j=1,..,n, is larger than -* max(underflow, 1/overflow). -* -* The bound on x(j) is also used to determine when a step in the -* columnwise method can be performed without fear of overflow. If -* the computed bound is greater than a large constant, x is scaled to -* prevent overflow, but if the bound overflows, x is set to 0, x(j) to -* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. -* -* Similarly, a row-wise scheme is used to solve A'*x = b. The basic -* algorithm for A upper triangular is -* -* for j = 1, ..., n -* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) -* end -* -* We simultaneously compute two bounds -* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j -* M(j) = bound on x(i), 1<=i<=j -* -* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we -* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. -* Then the bound on x(j) is -* -* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | -* -* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) -* 1<=i<=j -* -* and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater -* than max(underflow, 1/overflow). -* * ===================================================================== * * .. Parameters .. @@ -346,7 +431,7 @@ * ELSE * -* Compute the growth in A' * x = b. +* Compute the growth in A**T * x = b. * IF( UPPER ) THEN JFIRST = 1 @@ -554,7 +639,7 @@ * ELSE * -* Solve A' * x = b +* Solve A**T * x = b * DO 160 J = JFIRST, JLAST, JINC * @@ -666,7 +751,7 @@ ELSE * * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and -* scale = 0, and compute a solution to A'*x = 0. +* scale = 0, and compute a solution to A**T*x = 0. * DO 140 I = 1, N X( I ) = ZERO diff --git a/lib/linalg/drscl.f b/lib/linalg/drscl.f index 5a9ed00812..21ba19c11a 100644 --- a/lib/linalg/drscl.f +++ b/lib/linalg/drscl.f @@ -1,9 +1,93 @@ +*> \brief \b DRSCL multiplies a vector by the reciprocal of a real scalar. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DRSCL + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DRSCL( N, SA, SX, INCX ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* DOUBLE PRECISION SA +* .. +* .. Array Arguments .. +* DOUBLE PRECISION SX( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DRSCL multiplies an n-element real vector x by the real scalar 1/a. +*> This is done without overflow or underflow as long as +*> the final result x/a does not overflow or underflow. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of components of the vector x. +*> \endverbatim +*> +*> \param[in] SA +*> \verbatim +*> SA is DOUBLE PRECISION +*> The scalar a which is used to divide each component of x. +*> SA must be >= 0, or the subroutine will divide by zero. +*> \endverbatim +*> +*> \param[in,out] SX +*> \verbatim +*> SX is DOUBLE PRECISION array, dimension +*> (1+(N-1)*abs(INCX)) +*> The n-element vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between successive values of the vector SX. +*> > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= n +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleOTHERauxiliary +* +* ===================================================================== SUBROUTINE DRSCL( N, SA, SX, INCX ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. INTEGER INCX, N @@ -13,31 +97,6 @@ DOUBLE PRECISION SX( * ) * .. * -* Purpose -* ======= -* -* DRSCL multiplies an n-element real vector x by the real scalar 1/a. -* This is done without overflow or underflow as long as -* the final result x/a does not overflow or underflow. -* -* Arguments -* ========= -* -* N (input) INTEGER -* The number of components of the vector x. -* -* SA (input) DOUBLE PRECISION -* The scalar a which is used to divide each component of x. -* SA must be >= 0, or the subroutine will divide by zero. -* -* SX (input/output) DOUBLE PRECISION array, dimension -* (1+(N-1)*abs(INCX)) -* The n-element vector x. -* -* INCX (input) INTEGER -* The increment between successive values of the vector SX. -* > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= n -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/dscal.f b/lib/linalg/dscal.f index 81425a9474..3337de8e63 100644 --- a/lib/linalg/dscal.f +++ b/lib/linalg/dscal.f @@ -1,4 +1,63 @@ +*> \brief \b DSCAL +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DSCAL(N,DA,DX,INCX) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION DA +* INTEGER INCX,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DSCAL scales a vector by a constant. +*> uses unrolled loops for increment equal to one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 3/93 to return if incx .le. 0. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DSCAL(N,DA,DX,INCX) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. DOUBLE PRECISION DA INTEGER INCX,N @@ -7,19 +66,6 @@ DOUBLE PRECISION DX(*) * .. * -* Purpose -* ======= -* -* DSCAL scales a vector by a constant. -* uses unrolled loops for increment equal to one. -* -* Further Details -* =============== -* -* jack dongarra, linpack, 3/11/78. -* modified 3/93 to return if incx .le. 0. -* modified 12/3/93, array(1) declarations changed to array(*) -* * ===================================================================== * * .. Local Scalars .. @@ -29,34 +75,36 @@ INTRINSIC MOD * .. IF (N.LE.0 .OR. INCX.LE.0) RETURN - IF (INCX.EQ.1) GO TO 20 -* -* code for increment not equal to 1 -* - NINCX = N*INCX - DO 10 I = 1,NINCX,INCX - DX(I) = DA*DX(I) - 10 CONTINUE - RETURN + IF (INCX.EQ.1) THEN * * code for increment equal to 1 * * * clean-up loop * - 20 M = MOD(N,5) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DX(I) = DA*DX(I) - 30 CONTINUE - IF (N.LT.5) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,5 - DX(I) = DA*DX(I) - DX(I+1) = DA*DX(I+1) - DX(I+2) = DA*DX(I+2) - DX(I+3) = DA*DX(I+3) - DX(I+4) = DA*DX(I+4) - 50 CONTINUE + M = MOD(N,5) + IF (M.NE.0) THEN + DO I = 1,M + DX(I) = DA*DX(I) + END DO + IF (N.LT.5) RETURN + END IF + MP1 = M + 1 + DO I = MP1,N,5 + DX(I) = DA*DX(I) + DX(I+1) = DA*DX(I+1) + DX(I+2) = DA*DX(I+2) + DX(I+3) = DA*DX(I+3) + DX(I+4) = DA*DX(I+4) + END DO + ELSE +* +* code for increment not equal to 1 +* + NINCX = N*INCX + DO I = 1,NINCX,INCX + DX(I) = DA*DX(I) + END DO + END IF RETURN END diff --git a/lib/linalg/dswap.f b/lib/linalg/dswap.f index 3595fd8c10..e567bd93ec 100644 --- a/lib/linalg/dswap.f +++ b/lib/linalg/dswap.f @@ -1,4 +1,61 @@ +*> \brief \b DSWAP +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) +* +* .. Scalar Arguments .. +* INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*),DY(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> interchanges two vectors. +*> uses unrolled loops for increments equal one. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. INTEGER INCX,INCY,N * .. @@ -6,18 +63,6 @@ DOUBLE PRECISION DX(*),DY(*) * .. * -* Purpose -* ======= -* -* interchanges two vectors. -* uses unrolled loops for increments equal one. -* -* Further Details -* =============== -* -* jack dongarra, linpack, 3/11/78. -* modified 12/3/93, array(1) declarations changed to array(*) -* * ===================================================================== * * .. Local Scalars .. @@ -28,48 +73,50 @@ INTRINSIC MOD * .. IF (N.LE.0) RETURN - IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 -* -* code for unequal increments or equal increments not equal -* to 1 -* - IX = 1 - IY = 1 - IF (INCX.LT.0) IX = (-N+1)*INCX + 1 - IF (INCY.LT.0) IY = (-N+1)*INCY + 1 - DO 10 I = 1,N - DTEMP = DX(IX) - DX(IX) = DY(IY) - DY(IY) = DTEMP - IX = IX + INCX - IY = IY + INCY - 10 CONTINUE - RETURN + IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN * * code for both increments equal to 1 * * * clean-up loop * - 20 M = MOD(N,3) - IF (M.EQ.0) GO TO 40 - DO 30 I = 1,M - DTEMP = DX(I) - DX(I) = DY(I) - DY(I) = DTEMP - 30 CONTINUE - IF (N.LT.3) RETURN - 40 MP1 = M + 1 - DO 50 I = MP1,N,3 - DTEMP = DX(I) - DX(I) = DY(I) - DY(I) = DTEMP - DTEMP = DX(I+1) - DX(I+1) = DY(I+1) - DY(I+1) = DTEMP - DTEMP = DX(I+2) - DX(I+2) = DY(I+2) - DY(I+2) = DTEMP - 50 CONTINUE + M = MOD(N,3) + IF (M.NE.0) THEN + DO I = 1,M + DTEMP = DX(I) + DX(I) = DY(I) + DY(I) = DTEMP + END DO + IF (N.LT.3) RETURN + END IF + MP1 = M + 1 + DO I = MP1,N,3 + DTEMP = DX(I) + DX(I) = DY(I) + DY(I) = DTEMP + DTEMP = DX(I+1) + DX(I+1) = DY(I+1) + DY(I+1) = DTEMP + DTEMP = DX(I+2) + DX(I+2) = DY(I+2) + DY(I+2) = DTEMP + END DO + ELSE +* +* code for unequal increments or equal increments not equal +* to 1 +* + IX = 1 + IY = 1 + IF (INCX.LT.0) IX = (-N+1)*INCX + 1 + IF (INCY.LT.0) IY = (-N+1)*INCY + 1 + DO I = 1,N + DTEMP = DX(IX) + DX(IX) = DY(IY) + DY(IY) = DTEMP + IX = IX + INCX + IY = IY + INCY + END DO + END IF RETURN END diff --git a/lib/linalg/dtrmm.f b/lib/linalg/dtrmm.f index 931748dbeb..cbd5ce7034 100644 --- a/lib/linalg/dtrmm.f +++ b/lib/linalg/dtrmm.f @@ -1,4 +1,187 @@ +*> \brief \b DTRMM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA +* INTEGER LDA,LDB,M,N +* CHARACTER DIAG,SIDE,TRANSA,UPLO +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),B(LDB,*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRMM performs one of the matrix-matrix operations +*> +*> B := alpha*op( A )*B, or B := alpha*B*op( A ), +*> +*> where alpha is a scalar, B is an m by n matrix, A is a unit, or +*> non-unit, upper or lower triangular matrix and op( A ) is one of +*> +*> op( A ) = A or op( A ) = A**T. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> On entry, SIDE specifies whether op( A ) multiplies B from +*> the left or right as follows: +*> +*> SIDE = 'L' or 'l' B := alpha*op( A )*B. +*> +*> SIDE = 'R' or 'r' B := alpha*B*op( A ). +*> \endverbatim +*> +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> On entry, UPLO specifies whether the matrix A is an upper or +*> lower triangular matrix as follows: +*> +*> UPLO = 'U' or 'u' A is an upper triangular matrix. +*> +*> UPLO = 'L' or 'l' A is a lower triangular matrix. +*> \endverbatim +*> +*> \param[in] TRANSA +*> \verbatim +*> TRANSA is CHARACTER*1 +*> On entry, TRANSA specifies the form of op( A ) to be used in +*> the matrix multiplication as follows: +*> +*> TRANSA = 'N' or 'n' op( A ) = A. +*> +*> TRANSA = 'T' or 't' op( A ) = A**T. +*> +*> TRANSA = 'C' or 'c' op( A ) = A**T. +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> On entry, DIAG specifies whether or not A is unit triangular +*> as follows: +*> +*> DIAG = 'U' or 'u' A is assumed to be unit triangular. +*> +*> DIAG = 'N' or 'n' A is not assumed to be unit +*> triangular. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of B. M must be at +*> least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of B. N must be +*> at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. When alpha is +*> zero then A is not referenced and B need not be set before +*> entry. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m +*> when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. +*> Before entry with UPLO = 'U' or 'u', the leading k by k +*> upper triangular part of the array A must contain the upper +*> triangular matrix and the strictly lower triangular part of +*> A is not referenced. +*> Before entry with UPLO = 'L' or 'l', the leading k by k +*> lower triangular part of the array A must contain the lower +*> triangular matrix and the strictly upper triangular part of +*> A is not referenced. +*> Note that when DIAG = 'U' or 'u', the diagonal elements of +*> A are not referenced either, but are assumed to be unity. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. When SIDE = 'L' or 'l' then +*> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +*> then LDA must be at least max( 1, n ). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is DOUBLE PRECISION array of DIMENSION ( LDB, n ). +*> Before entry, the leading m by n part of the array B must +*> contain the matrix B, and on exit is overwritten by the +*> transformed matrix. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> On entry, LDB specifies the first dimension of B as declared +*> in the calling (sub) program. LDB must be at least +*> max( 1, m ). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level3 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 3 Blas routine. +*> +*> -- Written on 8-February-1989. +*> Jack Dongarra, Argonne National Laboratory. +*> Iain Duff, AERE Harwell. +*> Jeremy Du Croz, Numerical Algorithms Group Ltd. +*> Sven Hammarling, Numerical Algorithms Group Ltd. +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +* +* -- Reference BLAS level3 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER LDA,LDB,M,N @@ -8,123 +191,6 @@ DOUBLE PRECISION A(LDA,*),B(LDB,*) * .. * -* Purpose -* ======= -* -* DTRMM performs one of the matrix-matrix operations -* -* B := alpha*op( A )*B, or B := alpha*B*op( A ), -* -* where alpha is a scalar, B is an m by n matrix, A is a unit, or -* non-unit, upper or lower triangular matrix and op( A ) is one of -* -* op( A ) = A or op( A ) = A'. -* -* Arguments -* ========== -* -* SIDE - CHARACTER*1. -* On entry, SIDE specifies whether op( A ) multiplies B from -* the left or right as follows: -* -* SIDE = 'L' or 'l' B := alpha*op( A )*B. -* -* SIDE = 'R' or 'r' B := alpha*B*op( A ). -* -* Unchanged on exit. -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix A is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANSA - CHARACTER*1. -* On entry, TRANSA specifies the form of op( A ) to be used in -* the matrix multiplication as follows: -* -* TRANSA = 'N' or 'n' op( A ) = A. -* -* TRANSA = 'T' or 't' op( A ) = A'. -* -* TRANSA = 'C' or 'c' op( A ) = A'. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit triangular -* as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* M - INTEGER. -* On entry, M specifies the number of rows of B. M must be at -* least zero. -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the number of columns of B. N must be -* at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. When alpha is -* zero then A is not referenced and B need not be set before -* entry. -* Unchanged on exit. -* -* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m -* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. -* Before entry with UPLO = 'U' or 'u', the leading k by k -* upper triangular part of the array A must contain the upper -* triangular matrix and the strictly lower triangular part of -* A is not referenced. -* Before entry with UPLO = 'L' or 'l', the leading k by k -* lower triangular part of the array A must contain the lower -* triangular matrix and the strictly upper triangular part of -* A is not referenced. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced either, but are assumed to be unity. -* Unchanged on exit. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. When SIDE = 'L' or 'l' then -* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' -* then LDA must be at least max( 1, n ). -* Unchanged on exit. -* -* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). -* Before entry, the leading m by n part of the array B must -* contain the matrix B, and on exit is overwritten by the -* transformed matrix. -* -* LDB - INTEGER. -* On entry, LDB specifies the first dimension of B as declared -* in the calling (sub) program. LDB must be at least -* max( 1, m ). -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 3 Blas routine. -* -* -- Written on 8-February-1989. -* Jack Dongarra, Argonne National Laboratory. -* Iain Duff, AERE Harwell. -* Jeremy Du Croz, Numerical Algorithms Group Ltd. -* Sven Hammarling, Numerical Algorithms Group Ltd. -* * ===================================================================== * * .. External Functions .. @@ -234,7 +300,7 @@ END IF ELSE * -* Form B := alpha*A'*B. +* Form B := alpha*A**T*B. * IF (UPPER) THEN DO 110 J = 1,N @@ -300,7 +366,7 @@ END IF ELSE * -* Form B := alpha*B*A'. +* Form B := alpha*B*A**T. * IF (UPPER) THEN DO 260 K = 1,N diff --git a/lib/linalg/dtrmv.f b/lib/linalg/dtrmv.f index b1ca93a452..71459fe7c8 100644 --- a/lib/linalg/dtrmv.f +++ b/lib/linalg/dtrmv.f @@ -1,4 +1,157 @@ +*> \brief \b DTRMV +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) +* +* .. Scalar Arguments .. +* INTEGER INCX,LDA,N +* CHARACTER DIAG,TRANS,UPLO +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),X(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRMV performs one of the matrix-vector operations +*> +*> x := A*x, or x := A**T*x, +*> +*> where x is an n element vector and A is an n by n unit, or non-unit, +*> upper or lower triangular matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> On entry, UPLO specifies whether the matrix is an upper or +*> lower triangular matrix as follows: +*> +*> UPLO = 'U' or 'u' A is an upper triangular matrix. +*> +*> UPLO = 'L' or 'l' A is a lower triangular matrix. +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> On entry, TRANS specifies the operation to be performed as +*> follows: +*> +*> TRANS = 'N' or 'n' x := A*x. +*> +*> TRANS = 'T' or 't' x := A**T*x. +*> +*> TRANS = 'C' or 'c' x := A**T*x. +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> On entry, DIAG specifies whether or not A is unit +*> triangular as follows: +*> +*> DIAG = 'U' or 'u' A is assumed to be unit triangular. +*> +*> DIAG = 'N' or 'n' A is not assumed to be unit +*> triangular. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the order of the matrix A. +*> N must be at least zero. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, n ). +*> Before entry with UPLO = 'U' or 'u', the leading n by n +*> upper triangular part of the array A must contain the upper +*> triangular matrix and the strictly lower triangular part of +*> A is not referenced. +*> Before entry with UPLO = 'L' or 'l', the leading n by n +*> lower triangular part of the array A must contain the lower +*> triangular matrix and the strictly upper triangular part of +*> A is not referenced. +*> Note that when DIAG = 'U' or 'u', the diagonal elements of +*> A are not referenced either, but are assumed to be unity. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. LDA must be at least +*> max( 1, n ). +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is DOUBLE PRECISION array of dimension at least +*> ( 1 + ( n - 1 )*abs( INCX ) ). +*> Before entry, the incremented array X must contain the n +*> element vector x. On exit, X is overwritten with the +*> tranformed vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> On entry, INCX specifies the increment for the elements of +*> X. INCX must not be zero. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level2 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 2 Blas routine. +*> The vector and matrix arguments are not referenced when N = 0, or M = 0 +*> +*> -- Written on 22-October-1986. +*> Jack Dongarra, Argonne National Lab. +*> Jeremy Du Croz, Nag Central Office. +*> Sven Hammarling, Nag Central Office. +*> Richard Hanson, Sandia National Labs. +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) +* +* -- Reference BLAS level2 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. INTEGER INCX,LDA,N CHARACTER DIAG,TRANS,UPLO @@ -7,98 +160,6 @@ DOUBLE PRECISION A(LDA,*),X(*) * .. * -* Purpose -* ======= -* -* DTRMV performs one of the matrix-vector operations -* -* x := A*x, or x := A'*x, -* -* where x is an n element vector and A is an n by n unit, or non-unit, -* upper or lower triangular matrix. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the operation to be performed as -* follows: -* -* TRANS = 'N' or 'n' x := A*x. -* -* TRANS = 'T' or 't' x := A'*x. -* -* TRANS = 'C' or 'c' x := A'*x. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). -* Before entry with UPLO = 'U' or 'u', the leading n by n -* upper triangular part of the array A must contain the upper -* triangular matrix and the strictly lower triangular part of -* A is not referenced. -* Before entry with UPLO = 'L' or 'l', the leading n by n -* lower triangular part of the array A must contain the lower -* triangular matrix and the strictly upper triangular part of -* A is not referenced. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced either, but are assumed to be unity. -* Unchanged on exit. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. LDA must be at least -* max( 1, n ). -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. On exit, X is overwritten with the -* tranformed vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* * ===================================================================== * * .. Parameters .. @@ -221,7 +282,7 @@ END IF ELSE * -* Form x := A'*x. +* Form x := A**T*x. * IF (LSAME(UPLO,'U')) THEN IF (INCX.EQ.1) THEN diff --git a/lib/linalg/dtrsm.f b/lib/linalg/dtrsm.f index 888111042d..065df9a153 100644 --- a/lib/linalg/dtrsm.f +++ b/lib/linalg/dtrsm.f @@ -1,4 +1,191 @@ +*> \brief \b DTRSM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +* +* .. Scalar Arguments .. +* DOUBLE PRECISION ALPHA +* INTEGER LDA,LDB,M,N +* CHARACTER DIAG,SIDE,TRANSA,UPLO +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),B(LDB,*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRSM solves one of the matrix equations +*> +*> op( A )*X = alpha*B, or X*op( A ) = alpha*B, +*> +*> where alpha is a scalar, X and B are m by n matrices, A is a unit, or +*> non-unit, upper or lower triangular matrix and op( A ) is one of +*> +*> op( A ) = A or op( A ) = A**T. +*> +*> The matrix X is overwritten on B. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> On entry, SIDE specifies whether op( A ) appears on the left +*> or right of X as follows: +*> +*> SIDE = 'L' or 'l' op( A )*X = alpha*B. +*> +*> SIDE = 'R' or 'r' X*op( A ) = alpha*B. +*> \endverbatim +*> +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> On entry, UPLO specifies whether the matrix A is an upper or +*> lower triangular matrix as follows: +*> +*> UPLO = 'U' or 'u' A is an upper triangular matrix. +*> +*> UPLO = 'L' or 'l' A is a lower triangular matrix. +*> \endverbatim +*> +*> \param[in] TRANSA +*> \verbatim +*> TRANSA is CHARACTER*1 +*> On entry, TRANSA specifies the form of op( A ) to be used in +*> the matrix multiplication as follows: +*> +*> TRANSA = 'N' or 'n' op( A ) = A. +*> +*> TRANSA = 'T' or 't' op( A ) = A**T. +*> +*> TRANSA = 'C' or 'c' op( A ) = A**T. +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> On entry, DIAG specifies whether or not A is unit triangular +*> as follows: +*> +*> DIAG = 'U' or 'u' A is assumed to be unit triangular. +*> +*> DIAG = 'N' or 'n' A is not assumed to be unit +*> triangular. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> On entry, M specifies the number of rows of B. M must be at +*> least zero. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the number of columns of B. N must be +*> at least zero. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION. +*> On entry, ALPHA specifies the scalar alpha. When alpha is +*> zero then A is not referenced and B need not be set before +*> entry. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, k ), +*> where k is m when SIDE = 'L' or 'l' +*> and k is n when SIDE = 'R' or 'r'. +*> Before entry with UPLO = 'U' or 'u', the leading k by k +*> upper triangular part of the array A must contain the upper +*> triangular matrix and the strictly lower triangular part of +*> A is not referenced. +*> Before entry with UPLO = 'L' or 'l', the leading k by k +*> lower triangular part of the array A must contain the lower +*> triangular matrix and the strictly upper triangular part of +*> A is not referenced. +*> Note that when DIAG = 'U' or 'u', the diagonal elements of +*> A are not referenced either, but are assumed to be unity. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. When SIDE = 'L' or 'l' then +*> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +*> then LDA must be at least max( 1, n ). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is DOUBLE PRECISION array of DIMENSION ( LDB, n ). +*> Before entry, the leading m by n part of the array B must +*> contain the right-hand side matrix B, and on exit is +*> overwritten by the solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> On entry, LDB specifies the first dimension of B as declared +*> in the calling (sub) program. LDB must be at least +*> max( 1, m ). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level3 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Level 3 Blas routine. +*> +*> +*> -- Written on 8-February-1989. +*> Jack Dongarra, Argonne National Laboratory. +*> Iain Duff, AERE Harwell. +*> Jeremy Du Croz, Numerical Algorithms Group Ltd. +*> Sven Hammarling, Numerical Algorithms Group Ltd. +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) +* +* -- Reference BLAS level3 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER LDA,LDB,M,N @@ -8,126 +195,6 @@ DOUBLE PRECISION A(LDA,*),B(LDB,*) * .. * -* Purpose -* ======= -* -* DTRSM solves one of the matrix equations -* -* op( A )*X = alpha*B, or X*op( A ) = alpha*B, -* -* where alpha is a scalar, X and B are m by n matrices, A is a unit, or -* non-unit, upper or lower triangular matrix and op( A ) is one of -* -* op( A ) = A or op( A ) = A'. -* -* The matrix X is overwritten on B. -* -* Arguments -* ========== -* -* SIDE - CHARACTER*1. -* On entry, SIDE specifies whether op( A ) appears on the left -* or right of X as follows: -* -* SIDE = 'L' or 'l' op( A )*X = alpha*B. -* -* SIDE = 'R' or 'r' X*op( A ) = alpha*B. -* -* Unchanged on exit. -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix A is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANSA - CHARACTER*1. -* On entry, TRANSA specifies the form of op( A ) to be used in -* the matrix multiplication as follows: -* -* TRANSA = 'N' or 'n' op( A ) = A. -* -* TRANSA = 'T' or 't' op( A ) = A'. -* -* TRANSA = 'C' or 'c' op( A ) = A'. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit triangular -* as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* M - INTEGER. -* On entry, M specifies the number of rows of B. M must be at -* least zero. -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the number of columns of B. N must be -* at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. When alpha is -* zero then A is not referenced and B need not be set before -* entry. -* Unchanged on exit. -* -* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m -* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. -* Before entry with UPLO = 'U' or 'u', the leading k by k -* upper triangular part of the array A must contain the upper -* triangular matrix and the strictly lower triangular part of -* A is not referenced. -* Before entry with UPLO = 'L' or 'l', the leading k by k -* lower triangular part of the array A must contain the lower -* triangular matrix and the strictly upper triangular part of -* A is not referenced. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced either, but are assumed to be unity. -* Unchanged on exit. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. When SIDE = 'L' or 'l' then -* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' -* then LDA must be at least max( 1, n ). -* Unchanged on exit. -* -* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). -* Before entry, the leading m by n part of the array B must -* contain the right-hand side matrix B, and on exit is -* overwritten by the solution matrix X. -* -* LDB - INTEGER. -* On entry, LDB specifies the first dimension of B as declared -* in the calling (sub) program. LDB must be at least -* max( 1, m ). -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 3 Blas routine. -* -* -* -- Written on 8-February-1989. -* Jack Dongarra, Argonne National Laboratory. -* Iain Duff, AERE Harwell. -* Jeremy Du Croz, Numerical Algorithms Group Ltd. -* Sven Hammarling, Numerical Algorithms Group Ltd. -* * ===================================================================== * * .. External Functions .. @@ -243,7 +310,7 @@ END IF ELSE * -* Form B := alpha*inv( A' )*B. +* Form B := alpha*inv( A**T )*B. * IF (UPPER) THEN DO 130 J = 1,N @@ -319,7 +386,7 @@ END IF ELSE * -* Form B := alpha*B*inv( A' ). +* Form B := alpha*B*inv( A**T ). * IF (UPPER) THEN DO 310 K = N,1,-1 diff --git a/lib/linalg/dtrsv.f b/lib/linalg/dtrsv.f index 847fd3917d..e54303a93a 100644 --- a/lib/linalg/dtrsv.f +++ b/lib/linalg/dtrsv.f @@ -1,4 +1,153 @@ +*> \brief \b DTRSV +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) +* +* .. Scalar Arguments .. +* INTEGER INCX,LDA,N +* CHARACTER DIAG,TRANS,UPLO +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A(LDA,*),X(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRSV solves one of the systems of equations +*> +*> A*x = b, or A**T*x = b, +*> +*> where b and x are n element vectors and A is an n by n unit, or +*> non-unit, upper or lower triangular matrix. +*> +*> No test for singularity or near-singularity is included in this +*> routine. Such tests must be performed before calling this routine. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> On entry, UPLO specifies whether the matrix is an upper or +*> lower triangular matrix as follows: +*> +*> UPLO = 'U' or 'u' A is an upper triangular matrix. +*> +*> UPLO = 'L' or 'l' A is a lower triangular matrix. +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> On entry, TRANS specifies the equations to be solved as +*> follows: +*> +*> TRANS = 'N' or 'n' A*x = b. +*> +*> TRANS = 'T' or 't' A**T*x = b. +*> +*> TRANS = 'C' or 'c' A**T*x = b. +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> On entry, DIAG specifies whether or not A is unit +*> triangular as follows: +*> +*> DIAG = 'U' or 'u' A is assumed to be unit triangular. +*> +*> DIAG = 'N' or 'n' A is not assumed to be unit +*> triangular. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> On entry, N specifies the order of the matrix A. +*> N must be at least zero. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array of DIMENSION ( LDA, n ). +*> Before entry with UPLO = 'U' or 'u', the leading n by n +*> upper triangular part of the array A must contain the upper +*> triangular matrix and the strictly lower triangular part of +*> A is not referenced. +*> Before entry with UPLO = 'L' or 'l', the leading n by n +*> lower triangular part of the array A must contain the lower +*> triangular matrix and the strictly upper triangular part of +*> A is not referenced. +*> Note that when DIAG = 'U' or 'u', the diagonal elements of +*> A are not referenced either, but are assumed to be unity. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> On entry, LDA specifies the first dimension of A as declared +*> in the calling (sub) program. LDA must be at least +*> max( 1, n ). +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is DOUBLE PRECISION array of dimension at least +*> ( 1 + ( n - 1 )*abs( INCX ) ). +*> Before entry, the incremented array X must contain the n +*> element right-hand side vector b. On exit, X is overwritten +*> with the solution vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> On entry, INCX specifies the increment for the elements of +*> X. INCX must not be zero. +*> +*> Level 2 Blas routine. +*> +*> -- Written on 22-October-1986. +*> Jack Dongarra, Argonne National Lab. +*> Jeremy Du Croz, Nag Central Office. +*> Sven Hammarling, Nag Central Office. +*> Richard Hanson, Sandia National Labs. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup double_blas_level1 +* +* ===================================================================== SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. INTEGER INCX,LDA,N CHARACTER DIAG,TRANS,UPLO @@ -7,99 +156,6 @@ DOUBLE PRECISION A(LDA,*),X(*) * .. * -* Purpose -* ======= -* -* DTRSV solves one of the systems of equations -* -* A*x = b, or A'*x = b, -* -* where b and x are n element vectors and A is an n by n unit, or -* non-unit, upper or lower triangular matrix. -* -* No test for singularity or near-singularity is included in this -* routine. Such tests must be performed before calling this routine. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the equations to be solved as -* follows: -* -* TRANS = 'N' or 'n' A*x = b. -* -* TRANS = 'T' or 't' A'*x = b. -* -* TRANS = 'C' or 'c' A'*x = b. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). -* Before entry with UPLO = 'U' or 'u', the leading n by n -* upper triangular part of the array A must contain the upper -* triangular matrix and the strictly lower triangular part of -* A is not referenced. -* Before entry with UPLO = 'L' or 'l', the leading n by n -* lower triangular part of the array A must contain the lower -* triangular matrix and the strictly upper triangular part of -* A is not referenced. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced either, but are assumed to be unity. -* Unchanged on exit. -* -* LDA - INTEGER. -* On entry, LDA specifies the first dimension of A as declared -* in the calling (sub) program. LDA must be at least -* max( 1, n ). -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element right-hand side vector b. On exit, X is overwritten -* with the solution vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* * ===================================================================== * * .. Parameters .. @@ -221,7 +277,7 @@ END IF ELSE * -* Form x := inv( A' )*x. +* Form x := inv( A**T )*x. * IF (LSAME(UPLO,'U')) THEN IF (INCX.EQ.1) THEN diff --git a/lib/linalg/dtrti2.f b/lib/linalg/dtrti2.f index 2091b693dd..edf1b5b003 100644 --- a/lib/linalg/dtrti2.f +++ b/lib/linalg/dtrti2.f @@ -1,9 +1,119 @@ +*> \brief \b DTRTI2 computes the inverse of a triangular matrix (unblocked algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DTRTI2 + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER DIAG, UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRTI2 computes the inverse of a real upper or lower triangular +*> matrix. +*> +*> This is the Level 2 BLAS version of the algorithm. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the matrix A is upper or lower triangular. +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> Specifies whether or not the matrix A is unit triangular. +*> = 'N': Non-unit triangular +*> = 'U': Unit triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the triangular matrix A. If UPLO = 'U', the +*> leading n by n upper triangular part of the array A contains +*> the upper triangular matrix, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading n by n lower triangular part of the array A contains +*> the lower triangular matrix, and the strictly upper +*> triangular part of A is not referenced. If DIAG = 'U', the +*> diagonal elements of A are also not referenced and are +*> assumed to be 1. +*> +*> On exit, the (triangular) inverse of the original matrix, in +*> the same storage format. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-th argument had an illegal value +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date September 2012 +* +*> \ingroup doubleOTHERcomputational +* +* ===================================================================== SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) * -* -- LAPACK routine (version 3.2) -- +* -- LAPACK computational routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* September 2012 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO @@ -13,51 +123,6 @@ DOUBLE PRECISION A( LDA, * ) * .. * -* Purpose -* ======= -* -* DTRTI2 computes the inverse of a real upper or lower triangular -* matrix. -* -* This is the Level 2 BLAS version of the algorithm. -* -* Arguments -* ========= -* -* UPLO (input) CHARACTER*1 -* Specifies whether the matrix A is upper or lower triangular. -* = 'U': Upper triangular -* = 'L': Lower triangular -* -* DIAG (input) CHARACTER*1 -* Specifies whether or not the matrix A is unit triangular. -* = 'N': Non-unit triangular -* = 'U': Unit triangular -* -* N (input) INTEGER -* The order of the matrix A. N >= 0. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) -* On entry, the triangular matrix A. If UPLO = 'U', the -* leading n by n upper triangular part of the array A contains -* the upper triangular matrix, and the strictly lower -* triangular part of A is not referenced. If UPLO = 'L', the -* leading n by n lower triangular part of the array A contains -* the lower triangular matrix, and the strictly upper -* triangular part of A is not referenced. If DIAG = 'U', the -* diagonal elements of A are also not referenced and are -* assumed to be 1. -* -* On exit, the (triangular) inverse of the original matrix, in -* the same storage format. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -k, the k-th argument had an illegal value -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/dtrtri.f b/lib/linalg/dtrtri.f index 6a04cdf1f7..5d27ca56af 100644 --- a/lib/linalg/dtrtri.f +++ b/lib/linalg/dtrtri.f @@ -1,9 +1,118 @@ +*> \brief \b DTRTRI +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DTRTRI + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER DIAG, UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRTRI computes the inverse of a real upper or lower triangular +*> matrix A. +*> +*> This is the Level 3 BLAS version of the algorithm. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': A is upper triangular; +*> = 'L': A is lower triangular. +*> \endverbatim +*> +*> \param[in] DIAG +*> \verbatim +*> DIAG is CHARACTER*1 +*> = 'N': A is non-unit triangular; +*> = 'U': A is unit triangular. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the triangular matrix A. If UPLO = 'U', the +*> leading N-by-N upper triangular part of the array A contains +*> the upper triangular matrix, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading N-by-N lower triangular part of the array A contains +*> the lower triangular matrix, and the strictly upper +*> triangular part of A is not referenced. If DIAG = 'U', the +*> diagonal elements of A are also not referenced and are +*> assumed to be 1. +*> On exit, the (triangular) inverse of the original matrix, in +*> the same storage format. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, A(i,i) is exactly zero. The triangular +*> matrix is singular and its inverse can not be computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleOTHERcomputational +* +* ===================================================================== SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO ) * -* -- LAPACK routine (version 3.2) -- +* -- LAPACK computational routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO @@ -13,50 +122,6 @@ DOUBLE PRECISION A( LDA, * ) * .. * -* Purpose -* ======= -* -* DTRTRI computes the inverse of a real upper or lower triangular -* matrix A. -* -* This is the Level 3 BLAS version of the algorithm. -* -* Arguments -* ========= -* -* UPLO (input) CHARACTER*1 -* = 'U': A is upper triangular; -* = 'L': A is lower triangular. -* -* DIAG (input) CHARACTER*1 -* = 'N': A is non-unit triangular; -* = 'U': A is unit triangular. -* -* N (input) INTEGER -* The order of the matrix A. N >= 0. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) -* On entry, the triangular matrix A. If UPLO = 'U', the -* leading N-by-N upper triangular part of the array A contains -* the upper triangular matrix, and the strictly lower -* triangular part of A is not referenced. If UPLO = 'L', the -* leading N-by-N lower triangular part of the array A contains -* the lower triangular matrix, and the strictly upper -* triangular part of A is not referenced. If DIAG = 'U', the -* diagonal elements of A are also not referenced and are -* assumed to be 1. -* On exit, the (triangular) inverse of the original matrix, in -* the same storage format. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* > 0: if INFO = i, A(i,i) is exactly zero. The triangular -* matrix is singular and its inverse can not be computed. -* * ===================================================================== * * .. Parameters .. diff --git a/lib/linalg/idamax.f b/lib/linalg/idamax.f index fa45c1d69c..4233fcc273 100644 --- a/lib/linalg/idamax.f +++ b/lib/linalg/idamax.f @@ -1,4 +1,61 @@ +*> \brief \b IDAMAX +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* INTEGER FUNCTION IDAMAX(N,DX,INCX) +* +* .. Scalar Arguments .. +* INTEGER INCX,N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION DX(*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> IDAMAX finds the index of element having max. absolute value. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup aux_blas +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> jack dongarra, linpack, 3/11/78. +*> modified 3/93 to return if incx .le. 0. +*> modified 12/3/93, array(1) declarations changed to array(*) +*> \endverbatim +*> +* ===================================================================== INTEGER FUNCTION IDAMAX(N,DX,INCX) +* +* -- Reference BLAS level1 routine (version 3.4.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* * .. Scalar Arguments .. INTEGER INCX,N * .. @@ -6,18 +63,6 @@ DOUBLE PRECISION DX(*) * .. * -* Purpose -* ======= -* -* IDAMAX finds the index of element having max. absolute value. -* -* Further Details -* =============== -* -* jack dongarra, linpack, 3/11/78. -* modified 3/93 to return if incx .le. 0. -* modified 12/3/93, array(1) declarations changed to array(*) -* * ===================================================================== * * .. Local Scalars .. @@ -31,28 +76,31 @@ IF (N.LT.1 .OR. INCX.LE.0) RETURN IDAMAX = 1 IF (N.EQ.1) RETURN - IF (INCX.EQ.1) GO TO 20 -* -* code for increment not equal to 1 -* - IX = 1 - DMAX = DABS(DX(1)) - IX = IX + INCX - DO 10 I = 2,N - IF (DABS(DX(IX)).LE.DMAX) GO TO 5 - IDAMAX = I - DMAX = DABS(DX(IX)) - 5 IX = IX + INCX - 10 CONTINUE - RETURN + IF (INCX.EQ.1) THEN * * code for increment equal to 1 * - 20 DMAX = DABS(DX(1)) - DO 30 I = 2,N - IF (DABS(DX(I)).LE.DMAX) GO TO 30 - IDAMAX = I - DMAX = DABS(DX(I)) - 30 CONTINUE + DMAX = DABS(DX(1)) + DO I = 2,N + IF (DABS(DX(I)).GT.DMAX) THEN + IDAMAX = I + DMAX = DABS(DX(I)) + END IF + END DO + ELSE +* +* code for increment not equal to 1 +* + IX = 1 + DMAX = DABS(DX(1)) + IX = IX + INCX + DO I = 2,N + IF (DABS(DX(IX)).GT.DMAX) THEN + IDAMAX = I + DMAX = DABS(DX(IX)) + END IF + IX = IX + INCX + END DO + END IF RETURN END diff --git a/lib/linalg/ieeeck.f b/lib/linalg/ieeeck.f index b2c1d4909f..132e436770 100644 --- a/lib/linalg/ieeeck.f +++ b/lib/linalg/ieeeck.f @@ -1,43 +1,98 @@ +*> \brief \b IEEECK +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download IEEECK + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* INTEGER FUNCTION IEEECK( ISPEC, ZERO, ONE ) +* +* .. Scalar Arguments .. +* INTEGER ISPEC +* REAL ONE, ZERO +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> IEEECK is called from the ILAENV to verify that Infinity and +*> possibly NaN arithmetic is safe (i.e. will not trap). +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] ISPEC +*> \verbatim +*> ISPEC is INTEGER +*> Specifies whether to test just for inifinity arithmetic +*> or whether to test for infinity and NaN arithmetic. +*> = 0: Verify infinity arithmetic only. +*> = 1: Verify infinity and NaN arithmetic. +*> \endverbatim +*> +*> \param[in] ZERO +*> \verbatim +*> ZERO is REAL +*> Must contain the value 0.0 +*> This is passed to prevent the compiler from optimizing +*> away this code. +*> \endverbatim +*> +*> \param[in] ONE +*> \verbatim +*> ONE is REAL +*> Must contain the value 1.0 +*> This is passed to prevent the compiler from optimizing +*> away this code. +*> +*> RETURN VALUE: INTEGER +*> = 0: Arithmetic failed to produce the correct answers +*> = 1: Arithmetic produced the correct answers +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== INTEGER FUNCTION IEEECK( ISPEC, ZERO, ONE ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. INTEGER ISPEC REAL ONE, ZERO * .. * -* Purpose -* ======= -* -* IEEECK is called from the ILAENV to verify that Infinity and -* possibly NaN arithmetic is safe (i.e. will not trap). -* -* Arguments -* ========= -* -* ISPEC (input) INTEGER -* Specifies whether to test just for inifinity arithmetic -* or whether to test for infinity and NaN arithmetic. -* = 0: Verify infinity arithmetic only. -* = 1: Verify infinity and NaN arithmetic. -* -* ZERO (input) REAL -* Must contain the value 0.0 -* This is passed to prevent the compiler from optimizing -* away this code. -* -* ONE (input) REAL -* Must contain the value 1.0 -* This is passed to prevent the compiler from optimizing -* away this code. -* -* RETURN VALUE: INTEGER -* = 0: Arithmetic failed to produce the correct answers -* = 1: Arithmetic produced the correct answers +* ===================================================================== * * .. Local Scalars .. REAL NAN1, NAN2, NAN3, NAN4, NAN5, NAN6, NEGINF, @@ -112,7 +167,7 @@ * NAN5 = NEGINF*NEGZRO * - NAN6 = NAN5*0.0 + NAN6 = NAN5*ZERO * IF( NAN1.EQ.NAN1 ) THEN IEEECK = 0 diff --git a/lib/linalg/ilaenv.f b/lib/linalg/ilaenv.f index 8ed8aed657..867464de35 100644 --- a/lib/linalg/ilaenv.f +++ b/lib/linalg/ilaenv.f @@ -1,108 +1,177 @@ +*> \brief \b ILAENV +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ILAENV + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) +* +* .. Scalar Arguments .. +* CHARACTER*( * ) NAME, OPTS +* INTEGER ISPEC, N1, N2, N3, N4 +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ILAENV is called from the LAPACK routines to choose problem-dependent +*> parameters for the local environment. See ISPEC for a description of +*> the parameters. +*> +*> ILAENV returns an INTEGER +*> if ILAENV >= 0: ILAENV returns the value of the parameter specified by ISPEC +*> if ILAENV < 0: if ILAENV = -k, the k-th argument had an illegal value. +*> +*> This version provides a set of parameters which should give good, +*> but not optimal, performance on many of the currently available +*> computers. Users are encouraged to modify this subroutine to set +*> the tuning parameters for their particular machine using the option +*> and problem size information in the arguments. +*> +*> This routine will not function correctly if it is converted to all +*> lower case. Converting it to all upper case is allowed. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] ISPEC +*> \verbatim +*> ISPEC is INTEGER +*> Specifies the parameter to be returned as the value of +*> ILAENV. +*> = 1: the optimal blocksize; if this value is 1, an unblocked +*> algorithm will give the best performance. +*> = 2: the minimum block size for which the block routine +*> should be used; if the usable block size is less than +*> this value, an unblocked routine should be used. +*> = 3: the crossover point (in a block routine, for N less +*> than this value, an unblocked routine should be used) +*> = 4: the number of shifts, used in the nonsymmetric +*> eigenvalue routines (DEPRECATED) +*> = 5: the minimum column dimension for blocking to be used; +*> rectangular blocks must have dimension at least k by m, +*> where k is given by ILAENV(2,...) and m by ILAENV(5,...) +*> = 6: the crossover point for the SVD (when reducing an m by n +*> matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds +*> this value, a QR factorization is used first to reduce +*> the matrix to a triangular form.) +*> = 7: the number of processors +*> = 8: the crossover point for the multishift QR method +*> for nonsymmetric eigenvalue problems (DEPRECATED) +*> = 9: maximum size of the subproblems at the bottom of the +*> computation tree in the divide-and-conquer algorithm +*> (used by xGELSD and xGESDD) +*> =10: ieee NaN arithmetic can be trusted not to trap +*> =11: infinity arithmetic can be trusted not to trap +*> 12 <= ISPEC <= 16: +*> xHSEQR or one of its subroutines, +*> see IPARMQ for detailed explanation +*> \endverbatim +*> +*> \param[in] NAME +*> \verbatim +*> NAME is CHARACTER*(*) +*> The name of the calling subroutine, in either upper case or +*> lower case. +*> \endverbatim +*> +*> \param[in] OPTS +*> \verbatim +*> OPTS is CHARACTER*(*) +*> The character options to the subroutine NAME, concatenated +*> into a single character string. For example, UPLO = 'U', +*> TRANS = 'T', and DIAG = 'N' for a triangular routine would +*> be specified as OPTS = 'UTN'. +*> \endverbatim +*> +*> \param[in] N1 +*> \verbatim +*> N1 is INTEGER +*> \endverbatim +*> +*> \param[in] N2 +*> \verbatim +*> N2 is INTEGER +*> \endverbatim +*> +*> \param[in] N3 +*> \verbatim +*> N3 is INTEGER +*> \endverbatim +*> +*> \param[in] N4 +*> \verbatim +*> N4 is INTEGER +*> Problem dimensions for the subroutine NAME; these may not all +*> be required. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The following conventions have been used when calling ILAENV from the +*> LAPACK routines: +*> 1) OPTS is a concatenation of all of the character options to +*> subroutine NAME, in the same order that they appear in the +*> argument list for NAME, even if they are not used in determining +*> the value of the parameter specified by ISPEC. +*> 2) The problem dimensions N1, N2, N3, N4 are specified in the order +*> that they appear in the argument list for NAME. N1 is used +*> first, N2 second, and so on, and unused problem dimensions are +*> passed a value of -1. +*> 3) The parameter value returned by ILAENV is checked for validity in +*> the calling subroutine. For example, ILAENV is used to retrieve +*> the optimal blocksize for STRTRI as follows: +*> +*> NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) +*> IF( NB.LE.1 ) NB = MAX( 1, N ) +*> \endverbatim +*> +* ===================================================================== INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) * -* -- LAPACK auxiliary routine (version 3.2.1) -- -* -* -- April 2009 -- -* +* -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 * * .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER ISPEC, N1, N2, N3, N4 * .. * -* Purpose -* ======= -* -* ILAENV is called from the LAPACK routines to choose problem-dependent -* parameters for the local environment. See ISPEC for a description of -* the parameters. -* -* ILAENV returns an INTEGER -* if ILAENV >= 0: ILAENV returns the value of the parameter specified by ISPEC -* if ILAENV < 0: if ILAENV = -k, the k-th argument had an illegal value. -* -* This version provides a set of parameters which should give good, -* but not optimal, performance on many of the currently available -* computers. Users are encouraged to modify this subroutine to set -* the tuning parameters for their particular machine using the option -* and problem size information in the arguments. -* -* This routine will not function correctly if it is converted to all -* lower case. Converting it to all upper case is allowed. -* -* Arguments -* ========= -* -* ISPEC (input) INTEGER -* Specifies the parameter to be returned as the value of -* ILAENV. -* = 1: the optimal blocksize; if this value is 1, an unblocked -* algorithm will give the best performance. -* = 2: the minimum block size for which the block routine -* should be used; if the usable block size is less than -* this value, an unblocked routine should be used. -* = 3: the crossover point (in a block routine, for N less -* than this value, an unblocked routine should be used) -* = 4: the number of shifts, used in the nonsymmetric -* eigenvalue routines (DEPRECATED) -* = 5: the minimum column dimension for blocking to be used; -* rectangular blocks must have dimension at least k by m, -* where k is given by ILAENV(2,...) and m by ILAENV(5,...) -* = 6: the crossover point for the SVD (when reducing an m by n -* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds -* this value, a QR factorization is used first to reduce -* the matrix to a triangular form.) -* = 7: the number of processors -* = 8: the crossover point for the multishift QR method -* for nonsymmetric eigenvalue problems (DEPRECATED) -* = 9: maximum size of the subproblems at the bottom of the -* computation tree in the divide-and-conquer algorithm -* (used by xGELSD and xGESDD) -* =10: ieee NaN arithmetic can be trusted not to trap -* =11: infinity arithmetic can be trusted not to trap -* 12 <= ISPEC <= 16: -* xHSEQR or one of its subroutines, -* see IPARMQ for detailed explanation -* -* NAME (input) CHARACTER*(*) -* The name of the calling subroutine, in either upper case or -* lower case. -* -* OPTS (input) CHARACTER*(*) -* The character options to the subroutine NAME, concatenated -* into a single character string. For example, UPLO = 'U', -* TRANS = 'T', and DIAG = 'N' for a triangular routine would -* be specified as OPTS = 'UTN'. -* -* N1 (input) INTEGER -* N2 (input) INTEGER -* N3 (input) INTEGER -* N4 (input) INTEGER -* Problem dimensions for the subroutine NAME; these may not all -* be required. -* -* Further Details -* =============== -* -* The following conventions have been used when calling ILAENV from the -* LAPACK routines: -* 1) OPTS is a concatenation of all of the character options to -* subroutine NAME, in the same order that they appear in the -* argument list for NAME, even if they are not used in determining -* the value of the parameter specified by ISPEC. -* 2) The problem dimensions N1, N2, N3, N4 are specified in the order -* that they appear in the argument list for NAME. N1 is used -* first, N2 second, and so on, and unused problem dimensions are -* passed a value of -1. -* 3) The parameter value returned by ILAENV is checked for validity in -* the calling subroutine. For example, ILAENV is used to retrieve -* the optimal blocksize for STRTRI as follows: -* -* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) -* IF( NB.LE.1 ) NB = MAX( 1, N ) -* * ===================================================================== * * .. Local Scalars .. diff --git a/lib/linalg/iparmq.f b/lib/linalg/iparmq.f index 9a76099bb3..bd5bd7a0db 100644 --- a/lib/linalg/iparmq.f +++ b/lib/linalg/iparmq.f @@ -1,161 +1,229 @@ +*> \brief \b IPARMQ +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download IPARMQ + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* INTEGER FUNCTION IPARMQ( ISPEC, NAME, OPTS, N, ILO, IHI, LWORK ) +* +* .. Scalar Arguments .. +* INTEGER IHI, ILO, ISPEC, LWORK, N +* CHARACTER NAME*( * ), OPTS*( * ) +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> This program sets problem and machine dependent parameters +*> useful for xHSEQR and its subroutines. It is called whenever +*> ILAENV is called with 12 <= ISPEC <= 16 +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] ISPEC +*> \verbatim +*> ISPEC is integer scalar +*> ISPEC specifies which tunable parameter IPARMQ should +*> return. +*> +*> ISPEC=12: (INMIN) Matrices of order nmin or less +*> are sent directly to xLAHQR, the implicit +*> double shift QR algorithm. NMIN must be +*> at least 11. +*> +*> ISPEC=13: (INWIN) Size of the deflation window. +*> This is best set greater than or equal to +*> the number of simultaneous shifts NS. +*> Larger matrices benefit from larger deflation +*> windows. +*> +*> ISPEC=14: (INIBL) Determines when to stop nibbling and +*> invest in an (expensive) multi-shift QR sweep. +*> If the aggressive early deflation subroutine +*> finds LD converged eigenvalues from an order +*> NW deflation window and LD.GT.(NW*NIBBLE)/100, +*> then the next QR sweep is skipped and early +*> deflation is applied immediately to the +*> remaining active diagonal block. Setting +*> IPARMQ(ISPEC=14) = 0 causes TTQRE to skip a +*> multi-shift QR sweep whenever early deflation +*> finds a converged eigenvalue. Setting +*> IPARMQ(ISPEC=14) greater than or equal to 100 +*> prevents TTQRE from skipping a multi-shift +*> QR sweep. +*> +*> ISPEC=15: (NSHFTS) The number of simultaneous shifts in +*> a multi-shift QR iteration. +*> +*> ISPEC=16: (IACC22) IPARMQ is set to 0, 1 or 2 with the +*> following meanings. +*> 0: During the multi-shift QR sweep, +*> xLAQR5 does not accumulate reflections and +*> does not use matrix-matrix multiply to +*> update the far-from-diagonal matrix +*> entries. +*> 1: During the multi-shift QR sweep, +*> xLAQR5 and/or xLAQRaccumulates reflections and uses +*> matrix-matrix multiply to update the +*> far-from-diagonal matrix entries. +*> 2: During the multi-shift QR sweep. +*> xLAQR5 accumulates reflections and takes +*> advantage of 2-by-2 block structure during +*> matrix-matrix multiplies. +*> (If xTRMM is slower than xGEMM, then +*> IPARMQ(ISPEC=16)=1 may be more efficient than +*> IPARMQ(ISPEC=16)=2 despite the greater level of +*> arithmetic work implied by the latter choice.) +*> \endverbatim +*> +*> \param[in] NAME +*> \verbatim +*> NAME is character string +*> Name of the calling subroutine +*> \endverbatim +*> +*> \param[in] OPTS +*> \verbatim +*> OPTS is character string +*> This is a concatenation of the string arguments to +*> TTQRE. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is integer scalar +*> N is the order of the Hessenberg matrix H. +*> \endverbatim +*> +*> \param[in] ILO +*> \verbatim +*> ILO is INTEGER +*> \endverbatim +*> +*> \param[in] IHI +*> \verbatim +*> IHI is INTEGER +*> It is assumed that H is already upper triangular +*> in rows and columns 1:ILO-1 and IHI+1:N. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is integer scalar +*> The amount of workspace available. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Little is known about how best to choose these parameters. +*> It is possible to use different values of the parameters +*> for each of CHSEQR, DHSEQR, SHSEQR and ZHSEQR. +*> +*> It is probably best to choose different parameters for +*> different matrices and different parameters at different +*> times during the iteration, but this has not been +*> implemented --- yet. +*> +*> +*> The best choices of most of the parameters depend +*> in an ill-understood way on the relative execution +*> rate of xLAQR3 and xLAQR5 and on the nature of each +*> particular eigenvalue problem. Experiment may be the +*> only practical way to determine which choices are most +*> effective. +*> +*> Following is a list of default values supplied by IPARMQ. +*> These defaults may be adjusted in order to attain better +*> performance in any particular computational environment. +*> +*> IPARMQ(ISPEC=12) The xLAHQR vs xLAQR0 crossover point. +*> Default: 75. (Must be at least 11.) +*> +*> IPARMQ(ISPEC=13) Recommended deflation window size. +*> This depends on ILO, IHI and NS, the +*> number of simultaneous shifts returned +*> by IPARMQ(ISPEC=15). The default for +*> (IHI-ILO+1).LE.500 is NS. The default +*> for (IHI-ILO+1).GT.500 is 3*NS/2. +*> +*> IPARMQ(ISPEC=14) Nibble crossover point. Default: 14. +*> +*> IPARMQ(ISPEC=15) Number of simultaneous shifts, NS. +*> a multi-shift QR iteration. +*> +*> If IHI-ILO+1 is ... +*> +*> greater than ...but less ... the +*> or equal to ... than default is +*> +*> 0 30 NS = 2+ +*> 30 60 NS = 4+ +*> 60 150 NS = 10 +*> 150 590 NS = ** +*> 590 3000 NS = 64 +*> 3000 6000 NS = 128 +*> 6000 infinity NS = 256 +*> +*> (+) By default matrices of this order are +*> passed to the implicit double shift routine +*> xLAHQR. See IPARMQ(ISPEC=12) above. These +*> values of NS are used only in case of a rare +*> xLAHQR failure. +*> +*> (**) The asterisks (**) indicate an ad-hoc +*> function increasing from 10 to 64. +*> +*> IPARMQ(ISPEC=16) Select structured matrix multiply. +*> (See ISPEC=16 above for details.) +*> Default: 3. +*> \endverbatim +*> +* ===================================================================== INTEGER FUNCTION IPARMQ( ISPEC, NAME, OPTS, N, ILO, IHI, LWORK ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 -* +* November 2011 +* * .. Scalar Arguments .. INTEGER IHI, ILO, ISPEC, LWORK, N CHARACTER NAME*( * ), OPTS*( * ) * -* Purpose -* ======= -* -* This program sets problem and machine dependent parameters -* useful for xHSEQR and its subroutines. It is called whenever -* ILAENV is called with 12 <= ISPEC <= 16 -* -* Arguments -* ========= -* -* ISPEC (input) integer scalar -* ISPEC specifies which tunable parameter IPARMQ should -* return. -* -* ISPEC=12: (INMIN) Matrices of order nmin or less -* are sent directly to xLAHQR, the implicit -* double shift QR algorithm. NMIN must be -* at least 11. -* -* ISPEC=13: (INWIN) Size of the deflation window. -* This is best set greater than or equal to -* the number of simultaneous shifts NS. -* Larger matrices benefit from larger deflation -* windows. -* -* ISPEC=14: (INIBL) Determines when to stop nibbling and -* invest in an (expensive) multi-shift QR sweep. -* If the aggressive early deflation subroutine -* finds LD converged eigenvalues from an order -* NW deflation window and LD.GT.(NW*NIBBLE)/100, -* then the next QR sweep is skipped and early -* deflation is applied immediately to the -* remaining active diagonal block. Setting -* IPARMQ(ISPEC=14) = 0 causes TTQRE to skip a -* multi-shift QR sweep whenever early deflation -* finds a converged eigenvalue. Setting -* IPARMQ(ISPEC=14) greater than or equal to 100 -* prevents TTQRE from skipping a multi-shift -* QR sweep. -* -* ISPEC=15: (NSHFTS) The number of simultaneous shifts in -* a multi-shift QR iteration. -* -* ISPEC=16: (IACC22) IPARMQ is set to 0, 1 or 2 with the -* following meanings. -* 0: During the multi-shift QR sweep, -* xLAQR5 does not accumulate reflections and -* does not use matrix-matrix multiply to -* update the far-from-diagonal matrix -* entries. -* 1: During the multi-shift QR sweep, -* xLAQR5 and/or xLAQRaccumulates reflections and uses -* matrix-matrix multiply to update the -* far-from-diagonal matrix entries. -* 2: During the multi-shift QR sweep. -* xLAQR5 accumulates reflections and takes -* advantage of 2-by-2 block structure during -* matrix-matrix multiplies. -* (If xTRMM is slower than xGEMM, then -* IPARMQ(ISPEC=16)=1 may be more efficient than -* IPARMQ(ISPEC=16)=2 despite the greater level of -* arithmetic work implied by the latter choice.) -* -* NAME (input) character string -* Name of the calling subroutine -* -* OPTS (input) character string -* This is a concatenation of the string arguments to -* TTQRE. -* -* N (input) integer scalar -* N is the order of the Hessenberg matrix H. -* -* ILO (input) INTEGER -* IHI (input) INTEGER -* It is assumed that H is already upper triangular -* in rows and columns 1:ILO-1 and IHI+1:N. -* -* LWORK (input) integer scalar -* The amount of workspace available. -* -* Further Details -* =============== -* -* Little is known about how best to choose these parameters. -* It is possible to use different values of the parameters -* for each of CHSEQR, DHSEQR, SHSEQR and ZHSEQR. -* -* It is probably best to choose different parameters for -* different matrices and different parameters at different -* times during the iteration, but this has not been -* implemented --- yet. -* -* -* The best choices of most of the parameters depend -* in an ill-understood way on the relative execution -* rate of xLAQR3 and xLAQR5 and on the nature of each -* particular eigenvalue problem. Experiment may be the -* only practical way to determine which choices are most -* effective. -* -* Following is a list of default values supplied by IPARMQ. -* These defaults may be adjusted in order to attain better -* performance in any particular computational environment. -* -* IPARMQ(ISPEC=12) The xLAHQR vs xLAQR0 crossover point. -* Default: 75. (Must be at least 11.) -* -* IPARMQ(ISPEC=13) Recommended deflation window size. -* This depends on ILO, IHI and NS, the -* number of simultaneous shifts returned -* by IPARMQ(ISPEC=15). The default for -* (IHI-ILO+1).LE.500 is NS. The default -* for (IHI-ILO+1).GT.500 is 3*NS/2. -* -* IPARMQ(ISPEC=14) Nibble crossover point. Default: 14. -* -* IPARMQ(ISPEC=15) Number of simultaneous shifts, NS. -* a multi-shift QR iteration. -* -* If IHI-ILO+1 is ... -* -* greater than ...but less ... the -* or equal to ... than default is -* -* 0 30 NS = 2+ -* 30 60 NS = 4+ -* 60 150 NS = 10 -* 150 590 NS = ** -* 590 3000 NS = 64 -* 3000 6000 NS = 128 -* 6000 infinity NS = 256 -* -* (+) By default matrices of this order are -* passed to the implicit double shift routine -* xLAHQR. See IPARMQ(ISPEC=12) above. These -* values of NS are used only in case of a rare -* xLAHQR failure. -* -* (**) The asterisks (**) indicate an ad-hoc -* function increasing from 10 to 64. -* -* IPARMQ(ISPEC=16) Select structured matrix multiply. -* (See ISPEC=16 above for details.) -* Default: 3. -* -* ================================================================ +* ================================================================ * .. Parameters .. INTEGER INMIN, INWIN, INIBL, ISHFTS, IACC22 PARAMETER ( INMIN = 12, INWIN = 13, INIBL = 14, diff --git a/lib/linalg/lsame.f b/lib/linalg/lsame.f index 5e9047e14e..f19f9cda9e 100644 --- a/lib/linalg/lsame.f +++ b/lib/linalg/lsame.f @@ -1,83 +1,122 @@ - LOGICAL FUNCTION LSAME( CA, CB ) +*> \brief \b LSAME * -* -- LAPACK auxiliary routine (version 3.2) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* LOGICAL FUNCTION LSAME(CA,CB) +* +* .. Scalar Arguments .. +* CHARACTER CA,CB +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> LSAME returns .TRUE. if CA is the same letter as CB regardless of +*> case. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] CA +*> \verbatim +*> CA is CHARACTER*1 +*> \endverbatim +*> +*> \param[in] CB +*> \verbatim +*> CB is CHARACTER*1 +*> CA and CB specify the single characters to be compared. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup aux_blas +* +* ===================================================================== + LOGICAL FUNCTION LSAME(CA,CB) +* +* -- Reference BLAS level1 routine (version 3.1) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 * * .. Scalar Arguments .. - CHARACTER CA, CB + CHARACTER CA,CB * .. * -* Purpose -* ======= -* -* LSAME returns .TRUE. if CA is the same letter as CB regardless of -* case. -* -* Arguments -* ========= -* -* CA (input) CHARACTER*1 -* CB (input) CHARACTER*1 -* CA and CB specify the single characters to be compared. -* * ===================================================================== * * .. Intrinsic Functions .. - INTRINSIC ICHAR + INTRINSIC ICHAR * .. * .. Local Scalars .. - INTEGER INTA, INTB, ZCODE + INTEGER INTA,INTB,ZCODE * .. -* .. Executable Statements .. * * Test if the characters are equal * - LSAME = CA.EQ.CB - IF( LSAME ) - $ RETURN + LSAME = CA .EQ. CB + IF (LSAME) RETURN * * Now test for equivalence if both characters are alphabetic. * - ZCODE = ICHAR( 'Z' ) + ZCODE = ICHAR('Z') * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * - INTA = ICHAR( CA ) - INTB = ICHAR( CB ) + INTA = ICHAR(CA) + INTB = ICHAR(CB) * - IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN + IF (ZCODE.EQ.90 .OR. ZCODE.EQ.122) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * - IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 - IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 + IF (INTA.GE.97 .AND. INTA.LE.122) INTA = INTA - 32 + IF (INTB.GE.97 .AND. INTB.LE.122) INTB = INTB - 32 * - ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN + ELSE IF (ZCODE.EQ.233 .OR. ZCODE.EQ.169) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * - IF( INTA.GE.129 .AND. INTA.LE.137 .OR. - $ INTA.GE.145 .AND. INTA.LE.153 .OR. - $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 - IF( INTB.GE.129 .AND. INTB.LE.137 .OR. - $ INTB.GE.145 .AND. INTB.LE.153 .OR. - $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 + IF (INTA.GE.129 .AND. INTA.LE.137 .OR. + + INTA.GE.145 .AND. INTA.LE.153 .OR. + + INTA.GE.162 .AND. INTA.LE.169) INTA = INTA + 64 + IF (INTB.GE.129 .AND. INTB.LE.137 .OR. + + INTB.GE.145 .AND. INTB.LE.153 .OR. + + INTB.GE.162 .AND. INTB.LE.169) INTB = INTB + 64 * - ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN + ELSE IF (ZCODE.EQ.218 .OR. ZCODE.EQ.250) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * - IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 - IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 + IF (INTA.GE.225 .AND. INTA.LE.250) INTA = INTA - 32 + IF (INTB.GE.225 .AND. INTB.LE.250) INTB = INTB - 32 END IF - LSAME = INTA.EQ.INTB + LSAME = INTA .EQ. INTB * * RETURN * diff --git a/lib/linalg/xerbla.f b/lib/linalg/xerbla.f index 3a84150e32..3e93bc4e0e 100644 --- a/lib/linalg/xerbla.f +++ b/lib/linalg/xerbla.f @@ -1,34 +1,85 @@ +*> \brief \b XERBLA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download XERBLA + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE XERBLA( SRNAME, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER*(*) SRNAME +* INTEGER INFO +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> XERBLA is an error handler for the LAPACK routines. +*> It is called by an LAPACK routine if an input parameter has an +*> invalid value. A message is printed and execution stops. +*> +*> Installers may consider modifying the STOP statement in order to +*> call system-specific exception-handling facilities. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SRNAME +*> \verbatim +*> SRNAME is CHARACTER*(*) +*> The name of the routine which called XERBLA. +*> \endverbatim +*> +*> \param[in] INFO +*> \verbatim +*> INFO is INTEGER +*> The position of the invalid parameter in the parameter list +*> of the calling routine. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== SUBROUTINE XERBLA( SRNAME, INFO ) * -* -- LAPACK auxiliary routine (preliminary version) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 * * .. Scalar Arguments .. CHARACTER*(*) SRNAME INTEGER INFO * .. * -* Purpose -* ======= -* -* XERBLA is an error handler for the LAPACK routines. -* It is called by an LAPACK routine if an input parameter has an -* invalid value. A message is printed and execution stops. -* -* Installers may consider modifying the STOP statement in order to -* call system-specific exception-handling facilities. -* -* Arguments -* ========= -* -* SRNAME (input) CHARACTER*(*) -* The name of the routine which called XERBLA. -* -* INFO (input) INTEGER -* The position of the invalid parameter in the parameter list -* of the calling routine. -* * ===================================================================== * * .. Intrinsic Functions ..