ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9989 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2013-05-31 15:35:54 +00:00
parent b1430e190f
commit 2a70442035
35 changed files with 4008 additions and 2059 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
lib/linalg/Makefile.gfortran
View File

@ -1,23 +1,19 @@
# * # *
# *_________________________________________________________________________* # *_________________________________________________________________________*
# * Minimal BLAS/LAPACK Library for ATC # * Minimal BLAS/LAPACK Library for use by other LAMMPS packages
# 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
SHELL = /bin/sh SHELL = /bin/sh
# ------ FILES ------ # ------ 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.* README
FILES = $(SRC) Makefile.*
# ------ DEFINITIONS ------ # ------ DEFINITIONS ------
@ -29,7 +25,7 @@ OBJ = $(SRC:.f=.o)
FC = gfortran FC = gfortran
FFLAGS = -O3 -fPIC -march=native -mpc64 \ FFLAGS = -O3 -fPIC -march=native -mpc64 \
-ffast-math -funroll-loops -fstrict-aliasing -Wall -W -Wno-uninitialized -fno-second-underscore -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 -Wall -W -Wno-uninitialized -fno-second-underscore
ARCHIVE = ar ARCHIVE = ar
AR = ar AR = ar

18
lib/linalg/Makefile.mingw_cross
View File

@ -1,21 +1,17 @@
# * # *
# *_________________________________________________________________________* # *_________________________________________________________________________*
# * Minimal BLAS/LAPACK Library for ATC # * Minimal BLAS/LAPACK Library for use by other LAMMPS packages
# 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
SHELL = /bin/sh SHELL = /bin/sh
# ------ FILES ------ # ------ 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

10
lib/linalg/README
View File

@ -1,9 +1,13 @@
This directory has BLAS and LAPACK files needed by the USER-ATC and This directory has BLAS and LAPACK files needed by the USER-ATC and
USER-AWPMD packages, and possibly by other packages in the future. 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 Note that this is an *incomplete* subset of full BLAS/LAPACK.
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. 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 Build the library using one of the provided Makefile.* files or create
your own, specific to your compiler and system. For example: your own, specific to your compiler and system. For example:

119
lib/linalg/dasum.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
INTEGER INCX,N INTEGER INCX,N
* .. * ..
@ -6,18 +63,6 @@
DOUBLE PRECISION DX(*) 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 .. * .. Local Scalars ..
@ -30,33 +75,37 @@
DASUM = 0.0d0 DASUM = 0.0d0
DTEMP = 0.0d0 DTEMP = 0.0d0
IF (N.LE.0 .OR. INCX.LE.0) RETURN IF (N.LE.0 .OR. INCX.LE.0) RETURN
IF (INCX.EQ.1) GO TO 20 IF (INCX.EQ.1) THEN
*
* 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
*
* code for increment equal to 1 * code for increment equal to 1
* *
* *
* clean-up loop * clean-up loop
* *
20 M = MOD(N,6) M = MOD(N,6)
IF (M.EQ.0) GO TO 40 IF (M.NE.0) THEN
DO 30 I = 1,M DO I = 1,M
DTEMP = DTEMP + DABS(DX(I)) DTEMP = DTEMP + DABS(DX(I))
30 CONTINUE END DO
IF (N.LT.6) GO TO 60 IF (N.LT.6) THEN
40 MP1 = M + 1 DASUM = DTEMP
DO 50 I = MP1,N,6 RETURN
DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I+1)) + DABS(DX(I+2)) + END IF
+ DABS(DX(I+3)) + DABS(DX(I+4)) + DABS(DX(I+5)) END IF
50 CONTINUE MP1 = M + 1
60 DASUM = DTEMP 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 RETURN
END END

128
lib/linalg/daxpy.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
DOUBLE PRECISION DA DOUBLE PRECISION DA
INTEGER INCX,INCY,N INTEGER INCX,INCY,N
@ -7,18 +65,6 @@
DOUBLE PRECISION DX(*),DY(*) 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 .. * .. Local Scalars ..
@ -29,39 +75,41 @@
* .. * ..
IF (N.LE.0) RETURN IF (N.LE.0) RETURN
IF (DA.EQ.0.0d0) RETURN IF (DA.EQ.0.0d0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* 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
* *
* code for both increments equal to 1 * code for both increments equal to 1
* *
* *
* clean-up loop * clean-up loop
* *
20 M = MOD(N,4) M = MOD(N,4)
IF (M.EQ.0) GO TO 40 IF (M.NE.0) THEN
DO 30 I = 1,M DO I = 1,M
DY(I) = DY(I) + DA*DX(I) DY(I) = DY(I) + DA*DX(I)
30 CONTINUE END DO
IF (N.LT.4) RETURN END IF
40 MP1 = M + 1 IF (N.LT.4) RETURN
DO 50 I = MP1,N,4 MP1 = M + 1
DY(I) = DY(I) + DA*DX(I) DO I = MP1,N,4
DY(I+1) = DY(I+1) + DA*DX(I+1) DY(I) = DY(I) + DA*DX(I)
DY(I+2) = DY(I+2) + DA*DX(I+2) DY(I+1) = DY(I+1) + DA*DX(I+1)
DY(I+3) = DY(I+3) + DA*DX(I+3) DY(I+2) = DY(I+2) + DA*DX(I+2)
50 CONTINUE 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 RETURN
END END

133
lib/linalg/dcopy.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
INTEGER INCX,INCY,N INTEGER INCX,INCY,N
* .. * ..
@ -6,18 +63,6 @@
DOUBLE PRECISION DX(*),DY(*) 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 .. * .. Local Scalars ..
@ -27,42 +72,44 @@
INTRINSIC MOD INTRINSIC MOD
* .. * ..
IF (N.LE.0) RETURN IF (N.LE.0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* 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
* *
* code for both increments equal to 1 * code for both increments equal to 1
* *
* *
* clean-up loop * clean-up loop
* *
20 M = MOD(N,7) M = MOD(N,7)
IF (M.EQ.0) GO TO 40 IF (M.NE.0) THEN
DO 30 I = 1,M DO I = 1,M
DY(I) = DX(I) DY(I) = DX(I)
30 CONTINUE END DO
IF (N.LT.7) RETURN IF (N.LT.7) RETURN
40 MP1 = M + 1 END IF
DO 50 I = MP1,N,7 MP1 = M + 1
DY(I) = DX(I) DO I = MP1,N,7
DY(I+1) = DX(I+1) DY(I) = DX(I)
DY(I+2) = DX(I+2) DY(I+1) = DX(I+1)
DY(I+3) = DX(I+3) DY(I+2) = DX(I+2)
DY(I+4) = DX(I+4) DY(I+3) = DX(I+3)
DY(I+5) = DX(I+5) DY(I+4) = DX(I+4)
DY(I+6) = DX(I+6) DY(I+5) = DX(I+5)
50 CONTINUE 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 RETURN
END END

127
lib/linalg/ddot.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
INTEGER INCX,INCY,N INTEGER INCX,INCY,N
* .. * ..
@ -6,18 +63,6 @@
DOUBLE PRECISION DX(*),DY(*) 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 .. * .. Local Scalars ..
@ -30,39 +75,43 @@
DDOT = 0.0d0 DDOT = 0.0d0
DTEMP = 0.0d0 DTEMP = 0.0d0
IF (N.LE.0) RETURN IF (N.LE.0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* 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
* *
* code for both increments equal to 1 * code for both increments equal to 1
* *
* *
* clean-up loop * clean-up loop
* *
20 M = MOD(N,5) M = MOD(N,5)
IF (M.EQ.0) GO TO 40 IF (M.NE.0) THEN
DO 30 I = 1,M DO I = 1,M
DTEMP = DTEMP + DX(I)*DY(I) DTEMP = DTEMP + DX(I)*DY(I)
30 CONTINUE END DO
IF (N.LT.5) GO TO 60 IF (N.LT.5) THEN
40 MP1 = M + 1 DDOT=DTEMP
DO 50 I = MP1,N,5 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) + 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) $ DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
50 CONTINUE END DO
60 DDOT = DTEMP 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 RETURN
END END

179
lib/linalg/dgecon.f
View File

@ -1,12 +1,133 @@
*> \brief \b DGECON
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DGECON + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgecon.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgecon.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgecon.f">
*> [TXT]</a>
*> \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, SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
$ INFO ) $ INFO )
* *
* -- LAPACK routine (version 3.2) -- * -- LAPACK computational routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * November 2011
*
* Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM CHARACTER NORM
@ -18,52 +139,6 @@
DOUBLE PRECISION A( LDA, * ), WORK( * ) 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 .. * .. Parameters ..
@ -149,12 +224,12 @@
$ A, LDA, WORK, SU, WORK( 3*N+1 ), INFO ) $ A, LDA, WORK, SU, WORK( 3*N+1 ), INFO )
ELSE ELSE
* *
* Multiply by inv(U'). * Multiply by inv(U**T).
* *
CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A, CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
$ LDA, WORK, SU, WORK( 3*N+1 ), INFO ) $ 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, CALL DLATRS( 'Lower', 'Transpose', 'Unit', NORMIN, N, A,
$ LDA, WORK, SL, WORK( 2*N+1 ), INFO ) $ LDA, WORK, SL, WORK( 2*N+1 ), INFO )

320
lib/linalg/dgemm.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA DOUBLE PRECISION ALPHA,BETA
INTEGER K,LDA,LDB,LDC,M,N INTEGER K,LDA,LDB,LDC,M,N
@ -8,127 +201,6 @@
DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) 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 .. * .. External Functions ..
@ -249,7 +321,7 @@
90 CONTINUE 90 CONTINUE
ELSE ELSE
* *
* Form C := alpha*A'*B + beta*C * Form C := alpha*A**T*B + beta*C
* *
DO 120 J = 1,N DO 120 J = 1,N
DO 110 I = 1,M DO 110 I = 1,M
@ -268,7 +340,7 @@
ELSE ELSE
IF (NOTA) THEN IF (NOTA) THEN
* *
* Form C := alpha*A*B' + beta*C * Form C := alpha*A*B**T + beta*C
* *
DO 170 J = 1,N DO 170 J = 1,N
IF (BETA.EQ.ZERO) THEN IF (BETA.EQ.ZERO) THEN
@ -291,7 +363,7 @@
170 CONTINUE 170 CONTINUE
ELSE ELSE
* *
* Form C := alpha*A'*B' + beta*C * Form C := alpha*A**T*B**T + beta*C
* *
DO 200 J = 1,N DO 200 J = 1,N
DO 190 I = 1,M DO 190 I = 1,M

256
lib/linalg/dgemv.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA DOUBLE PRECISION ALPHA,BETA
INTEGER INCX,INCY,LDA,M,N INTEGER INCX,INCY,LDA,M,N
@ -8,98 +170,6 @@
DOUBLE PRECISION A(LDA,*),X(*),Y(*) 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 .. * .. Parameters ..
@ -231,7 +301,7 @@
END IF END IF
ELSE ELSE
* *
* Form y := alpha*A'*x + y. * Form y := alpha*A**T*x + y.
* *
JY = KY JY = KY
IF (INCX.EQ.1) THEN IF (INCX.EQ.1) THEN

207
lib/linalg/dger.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
DOUBLE PRECISION ALPHA DOUBLE PRECISION ALPHA
INTEGER INCX,INCY,LDA,M,N INTEGER INCX,INCY,LDA,M,N
@ -7,77 +143,6 @@
DOUBLE PRECISION A(LDA,*),X(*),Y(*) 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 .. * .. Parameters ..

155
lib/linalg/dgetf2.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgetf2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgetf2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgetf2.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INFO, LDA, M, N INTEGER INFO, LDA, M, N
@ -13,49 +121,6 @@
DOUBLE PRECISION A( LDA, * ) 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 .. * .. Parameters ..

155
lib/linalg/dgetrf.f
View File

@ -1,9 +1,117 @@
*> \brief \b DGETRF
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DGETRF + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgetrf.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgetrf.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgetrf.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INFO, LDA, M, N INTEGER INFO, LDA, M, N
@ -13,49 +121,6 @@
DOUBLE PRECISION A( LDA, * ) 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 .. * .. Parameters ..

164
lib/linalg/dgetri.f
View File

@ -1,9 +1,123 @@
*> \brief \b DGETRI
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DGETRI + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgetri.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgetri.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgetri.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, N INTEGER INFO, LDA, LWORK, N
@ -13,52 +127,6 @@
DOUBLE PRECISION A( LDA, * ), WORK( * ) 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 .. * .. Parameters ..

103
lib/linalg/dlabad.f
View File

@ -1,39 +1,88 @@
*> \brief \b DLABAD
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLABAD + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlabad.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlabad.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlabad.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
DOUBLE PRECISION LARGE, SMALL 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 .. * .. Intrinsic Functions ..

197
lib/linalg/dlacn2.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlacn2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlacn2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlacn2.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER KASE, N INTEGER KASE, N
@ -14,63 +150,6 @@
DOUBLE PRECISION V( * ), X( * ) 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 .. * .. Parameters ..

184
lib/linalg/dlange.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlange.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlange.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlange.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM CHARACTER NORM
@ -13,56 +127,6 @@
DOUBLE PRECISION A( LDA, * ), WORK( * ) 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 .. * .. Parameters ..
@ -71,17 +135,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J INTEGER I, J
DOUBLE PRECISION SCALE, SUM, VALUE DOUBLE PRECISION SCALE, SUM, VALUE, TEMP
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL DLASSQ EXTERNAL DLASSQ
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME LOGICAL LSAME, DISNAN
EXTERNAL LSAME EXTERNAL LSAME, DISNAN
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT INTRINSIC ABS, MIN, SQRT
* .. * ..
* .. Executable Statements .. * .. Executable Statements ..
* *
@ -94,7 +158,8 @@
VALUE = ZERO VALUE = ZERO
DO 20 J = 1, N DO 20 J = 1, N
DO 10 I = 1, M 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 10 CONTINUE
20 CONTINUE 20 CONTINUE
ELSE IF( ( LSAME( NORM, 'O' ) ) .OR. ( NORM.EQ.'1' ) ) THEN ELSE IF( ( LSAME( NORM, 'O' ) ) .OR. ( NORM.EQ.'1' ) ) THEN
@ -107,7 +172,7 @@
DO 30 I = 1, M DO 30 I = 1, M
SUM = SUM + ABS( A( I, J ) ) SUM = SUM + ABS( A( I, J ) )
30 CONTINUE 30 CONTINUE
VALUE = MAX( VALUE, SUM ) IF( VALUE.LT.SUM .OR. DISNAN( SUM ) ) VALUE = SUM
40 CONTINUE 40 CONTINUE
ELSE IF( LSAME( NORM, 'I' ) ) THEN ELSE IF( LSAME( NORM, 'I' ) ) THEN
* *
@ -123,7 +188,8 @@
70 CONTINUE 70 CONTINUE
VALUE = ZERO VALUE = ZERO
DO 80 I = 1, M DO 80 I = 1, M
VALUE = MAX( VALUE, WORK( I ) ) TEMP = WORK( I )
IF( VALUE.LT.TEMP .OR. DISNAN( TEMP ) ) VALUE = TEMP
80 CONTINUE 80 CONTINUE
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *

156
lib/linalg/dlassq.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlassq.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlassq.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlassq.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INCX, N INTEGER INCX, N
@ -13,47 +116,6 @@
DOUBLE PRECISION X( * ) 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 .. * .. Parameters ..
@ -64,6 +126,10 @@
INTEGER IX INTEGER IX
DOUBLE PRECISION ABSXI DOUBLE PRECISION ABSXI
* .. * ..
* .. External Functions ..
LOGICAL DISNAN
EXTERNAL DISNAN
* ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS INTRINSIC ABS
* .. * ..
@ -71,8 +137,8 @@
* *
IF( N.GT.0 ) THEN IF( N.GT.0 ) THEN
DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX 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 IF( SCALE.LT.ABSXI ) THEN
SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2 SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2
SCALE = ABSXI SCALE = ABSXI

161
lib/linalg/dlaswp.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaswp.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaswp.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaswp.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INCX, K1, K2, LDA, N INTEGER INCX, K1, K2, LDA, N
@ -13,49 +127,6 @@
DOUBLE PRECISION A( LDA, * ) 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 .. * .. Local Scalars ..

399
lib/linalg/dlatrs.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlatrs.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlatrs.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlatrs.f">
*> [TXT]</a>
*> \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, SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
$ CNORM, INFO ) $ 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER DIAG, NORMIN, TRANS, UPLO CHARACTER DIAG, NORMIN, TRANS, UPLO
@ -15,158 +252,6 @@
DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * ) 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 .. * .. Parameters ..
@ -346,7 +431,7 @@
* *
ELSE ELSE
* *
* Compute the growth in A' * x = b. * Compute the growth in A**T * x = b.
* *
IF( UPPER ) THEN IF( UPPER ) THEN
JFIRST = 1 JFIRST = 1
@ -554,7 +639,7 @@
* *
ELSE ELSE
* *
* Solve A' * x = b * Solve A**T * x = b
* *
DO 160 J = JFIRST, JLAST, JINC DO 160 J = JFIRST, JLAST, JINC
* *
@ -666,7 +751,7 @@
ELSE ELSE
* *
* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and * 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 DO 140 I = 1, N
X( I ) = ZERO X( I ) = ZERO

113
lib/linalg/drscl.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/drscl.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/drscl.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/drscl.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INCX, N INTEGER INCX, N
@ -13,31 +97,6 @@
DOUBLE PRECISION SX( * ) 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 .. * .. Parameters ..

120
lib/linalg/dscal.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
DOUBLE PRECISION DA DOUBLE PRECISION DA
INTEGER INCX,N INTEGER INCX,N
@ -7,19 +66,6 @@
DOUBLE PRECISION DX(*) 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 .. * .. Local Scalars ..
@ -29,34 +75,36 @@
INTRINSIC MOD INTRINSIC MOD
* .. * ..
IF (N.LE.0 .OR. INCX.LE.0) RETURN IF (N.LE.0 .OR. INCX.LE.0) RETURN
IF (INCX.EQ.1) GO TO 20 IF (INCX.EQ.1) THEN
*
* code for increment not equal to 1
*
NINCX = N*INCX
DO 10 I = 1,NINCX,INCX
DX(I) = DA*DX(I)
10 CONTINUE
RETURN
* *
* code for increment equal to 1 * code for increment equal to 1
* *
* *
* clean-up loop * clean-up loop
* *
20 M = MOD(N,5) M = MOD(N,5)
IF (M.EQ.0) GO TO 40 IF (M.NE.0) THEN
DO 30 I = 1,M DO I = 1,M
DX(I) = DA*DX(I) DX(I) = DA*DX(I)
30 CONTINUE END DO
IF (N.LT.5) RETURN IF (N.LT.5) RETURN
40 MP1 = M + 1 END IF
DO 50 I = MP1,N,5 MP1 = M + 1
DX(I) = DA*DX(I) DO I = MP1,N,5
DX(I+1) = DA*DX(I+1) DX(I) = DA*DX(I)
DX(I+2) = DA*DX(I+2) DX(I+1) = DA*DX(I+1)
DX(I+3) = DA*DX(I+3) DX(I+2) = DA*DX(I+2)
DX(I+4) = DA*DX(I+4) DX(I+3) = DA*DX(I+3)
50 CONTINUE 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 RETURN
END END

145
lib/linalg/dswap.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
INTEGER INCX,INCY,N INTEGER INCX,INCY,N
* .. * ..
@ -6,18 +63,6 @@
DOUBLE PRECISION DX(*),DY(*) 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 .. * .. Local Scalars ..
@ -28,48 +73,50 @@
INTRINSIC MOD INTRINSIC MOD
* .. * ..
IF (N.LE.0) RETURN IF (N.LE.0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* code for unequal increments or equal increments not equal
* to 1
*
IX = 1
IY = 1
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DTEMP = DX(IX)
DX(IX) = DY(IY)
DY(IY) = DTEMP
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
* *
* code for both increments equal to 1 * code for both increments equal to 1
* *
* *
* clean-up loop * clean-up loop
* *
20 M = MOD(N,3) M = MOD(N,3)
IF (M.EQ.0) GO TO 40 IF (M.NE.0) THEN
DO 30 I = 1,M DO I = 1,M
DTEMP = DX(I) DTEMP = DX(I)
DX(I) = DY(I) DX(I) = DY(I)
DY(I) = DTEMP DY(I) = DTEMP
30 CONTINUE END DO
IF (N.LT.3) RETURN IF (N.LT.3) RETURN
40 MP1 = M + 1 END IF
DO 50 I = MP1,N,3 MP1 = M + 1
DTEMP = DX(I) DO I = MP1,N,3
DX(I) = DY(I) DTEMP = DX(I)
DY(I) = DTEMP DX(I) = DY(I)
DTEMP = DX(I+1) DY(I) = DTEMP
DX(I+1) = DY(I+1) DTEMP = DX(I+1)
DY(I+1) = DTEMP DX(I+1) = DY(I+1)
DTEMP = DX(I+2) DY(I+1) = DTEMP
DX(I+2) = DY(I+2) DTEMP = DX(I+2)
DY(I+2) = DTEMP DX(I+2) = DY(I+2)
50 CONTINUE 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 RETURN
END END

304
lib/linalg/dtrmm.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
DOUBLE PRECISION ALPHA DOUBLE PRECISION ALPHA
INTEGER LDA,LDB,M,N INTEGER LDA,LDB,M,N
@ -8,123 +191,6 @@
DOUBLE PRECISION A(LDA,*),B(LDB,*) 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 .. * .. External Functions ..
@ -234,7 +300,7 @@
END IF END IF
ELSE ELSE
* *
* Form B := alpha*A'*B. * Form B := alpha*A**T*B.
* *
IF (UPPER) THEN IF (UPPER) THEN
DO 110 J = 1,N DO 110 J = 1,N
@ -300,7 +366,7 @@
END IF END IF
ELSE ELSE
* *
* Form B := alpha*B*A'. * Form B := alpha*B*A**T.
* *
IF (UPPER) THEN IF (UPPER) THEN
DO 260 K = 1,N DO 260 K = 1,N

247
lib/linalg/dtrmv.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
INTEGER INCX,LDA,N INTEGER INCX,LDA,N
CHARACTER DIAG,TRANS,UPLO CHARACTER DIAG,TRANS,UPLO
@ -7,98 +160,6 @@
DOUBLE PRECISION A(LDA,*),X(*) 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 .. * .. Parameters ..
@ -221,7 +282,7 @@
END IF END IF
ELSE ELSE
* *
* Form x := A'*x. * Form x := A**T*x.
* *
IF (LSAME(UPLO,'U')) THEN IF (LSAME(UPLO,'U')) THEN
IF (INCX.EQ.1) THEN IF (INCX.EQ.1) THEN

311
lib/linalg/dtrsm.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
DOUBLE PRECISION ALPHA DOUBLE PRECISION ALPHA
INTEGER LDA,LDB,M,N INTEGER LDA,LDB,M,N
@ -8,126 +195,6 @@
DOUBLE PRECISION A(LDA,*),B(LDB,*) 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 .. * .. External Functions ..
@ -243,7 +310,7 @@
END IF END IF
ELSE ELSE
* *
* Form B := alpha*inv( A' )*B. * Form B := alpha*inv( A**T )*B.
* *
IF (UPPER) THEN IF (UPPER) THEN
DO 130 J = 1,N DO 130 J = 1,N
@ -319,7 +386,7 @@
END IF END IF
ELSE ELSE
* *
* Form B := alpha*B*inv( A' ). * Form B := alpha*B*inv( A**T ).
* *
IF (UPPER) THEN IF (UPPER) THEN
DO 310 K = N,1,-1 DO 310 K = N,1,-1

244
lib/linalg/dtrsv.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
INTEGER INCX,LDA,N INTEGER INCX,LDA,N
CHARACTER DIAG,TRANS,UPLO CHARACTER DIAG,TRANS,UPLO
@ -7,99 +156,6 @@
DOUBLE PRECISION A(LDA,*),X(*) 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 .. * .. Parameters ..
@ -221,7 +277,7 @@
END IF END IF
ELSE ELSE
* *
* Form x := inv( A' )*x. * Form x := inv( A**T )*x.
* *
IF (LSAME(UPLO,'U')) THEN IF (LSAME(UPLO,'U')) THEN
IF (INCX.EQ.1) THEN IF (INCX.EQ.1) THEN

159
lib/linalg/dtrti2.f
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrti2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrti2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrti2.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * September 2012
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER DIAG, UPLO CHARACTER DIAG, UPLO
@ -13,51 +123,6 @@
DOUBLE PRECISION A( LDA, * ) 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 .. * .. Parameters ..

157
lib/linalg/dtrtri.f
View File

@ -1,9 +1,118 @@
*> \brief \b DTRTRI
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DTRTRI + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrtri.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrtri.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrtri.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER DIAG, UPLO CHARACTER DIAG, UPLO
@ -13,50 +122,6 @@
DOUBLE PRECISION A( LDA, * ) 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 .. * .. Parameters ..

112
lib/linalg/idamax.f
View File

@ -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) 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 .. * .. Scalar Arguments ..
INTEGER INCX,N INTEGER INCX,N
* .. * ..
@ -6,18 +63,6 @@
DOUBLE PRECISION DX(*) 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 .. * .. Local Scalars ..
@ -31,28 +76,31 @@
IF (N.LT.1 .OR. INCX.LE.0) RETURN IF (N.LT.1 .OR. INCX.LE.0) RETURN
IDAMAX = 1 IDAMAX = 1
IF (N.EQ.1) RETURN IF (N.EQ.1) RETURN
IF (INCX.EQ.1) GO TO 20 IF (INCX.EQ.1) THEN
*
* 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
* *
* code for increment equal to 1 * code for increment equal to 1
* *
20 DMAX = DABS(DX(1)) DMAX = DABS(DX(1))
DO 30 I = 2,N DO I = 2,N
IF (DABS(DX(I)).LE.DMAX) GO TO 30 IF (DABS(DX(I)).GT.DMAX) THEN
IDAMAX = I IDAMAX = I
DMAX = DABS(DX(I)) DMAX = DABS(DX(I))
30 CONTINUE 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 RETURN
END END

117
lib/linalg/ieeeck.f
View File

@ -1,43 +1,98 @@
*> \brief \b IEEECK
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download IEEECK + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ieeeck.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ieeeck.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ieeeck.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER ISPEC INTEGER ISPEC
REAL ONE, ZERO 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 .. * .. Local Scalars ..
REAL NAN1, NAN2, NAN3, NAN4, NAN5, NAN6, NEGINF, REAL NAN1, NAN2, NAN3, NAN4, NAN5, NAN6, NEGINF,
@ -112,7 +167,7 @@
* *
NAN5 = NEGINF*NEGZRO NAN5 = NEGINF*NEGZRO
* *
NAN6 = NAN5*0.0 NAN6 = NAN5*ZERO
* *
IF( NAN1.EQ.NAN1 ) THEN IF( NAN1.EQ.NAN1 ) THEN
IEEECK = 0 IEEECK = 0

259
lib/linalg/ilaenv.f
View File

@ -1,108 +1,177 @@
*> \brief \b ILAENV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ILAENV + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ilaenv.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ilaenv.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ilaenv.f">
*> [TXT]</a>
*> \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 ) INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 )
* *
* -- LAPACK auxiliary routine (version 3.2.1) -- * -- LAPACK auxiliary routine (version 3.4.0) --
*
* -- April 2009 --
*
* -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER*( * ) NAME, OPTS CHARACTER*( * ) NAME, OPTS
INTEGER ISPEC, N1, N2, N3, N4 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 .. * .. Local Scalars ..

368
lib/linalg/iparmq.f
View File

@ -1,161 +1,229 @@
*> \brief \b IPARMQ
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download IPARMQ + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/iparmq.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/iparmq.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/iparmq.f">
*> [TXT]</a>
*> \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 ) 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, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006 * November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER IHI, ILO, ISPEC, LWORK, N INTEGER IHI, ILO, ISPEC, LWORK, N
CHARACTER NAME*( * ), OPTS*( * ) 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 .. * .. Parameters ..
INTEGER INMIN, INWIN, INIBL, ISHFTS, IACC22 INTEGER INMIN, INWIN, INIBL, ISHFTS, IACC22
PARAMETER ( INMIN = 12, INWIN = 13, INIBL = 14, PARAMETER ( INMIN = 12, INWIN = 13, INIBL = 14,

121
lib/linalg/lsame.f
View File

@ -1,83 +1,122 @@
LOGICAL FUNCTION LSAME( CA, CB ) *> \brief \b LSAME
* *
* -- LAPACK auxiliary routine (version 3.2) -- * =========== DOCUMENTATION ===========
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *
* November 2006 * 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 .. * .. 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 Functions ..
INTRINSIC ICHAR INTRINSIC ICHAR
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER INTA, INTB, ZCODE INTEGER INTA,INTB,ZCODE
* .. * ..
* .. Executable Statements ..
* *
* Test if the characters are equal * Test if the characters are equal
* *
LSAME = CA.EQ.CB LSAME = CA .EQ. CB
IF( LSAME ) IF (LSAME) RETURN
$ RETURN
* *
* Now test for equivalence if both characters are alphabetic. * 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 * Use 'Z' rather than 'A' so that ASCII can be detected on Prime
* machines, on which ICHAR returns a value with bit 8 set. * 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 Prime machines returns 193 which is the same as
* ICHAR('A') on an EBCDIC machine. * ICHAR('A') on an EBCDIC machine.
* *
INTA = ICHAR( CA ) INTA = ICHAR(CA)
INTB = ICHAR( CB ) 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 * ASCII is assumed - ZCODE is the ASCII code of either lower or
* upper case 'Z'. * upper case 'Z'.
* *
IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF (INTA.GE.97 .AND. INTA.LE.122) INTA = INTA - 32
IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 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 * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
* upper case 'Z'. * upper case 'Z'.
* *
IF( INTA.GE.129 .AND. INTA.LE.137 .OR. IF (INTA.GE.129 .AND. INTA.LE.137 .OR.
$ INTA.GE.145 .AND. INTA.LE.153 .OR. + INTA.GE.145 .AND. INTA.LE.153 .OR.
$ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 + INTA.GE.162 .AND. INTA.LE.169) INTA = INTA + 64
IF( INTB.GE.129 .AND. INTB.LE.137 .OR. IF (INTB.GE.129 .AND. INTB.LE.137 .OR.
$ INTB.GE.145 .AND. INTB.LE.153 .OR. + INTB.GE.145 .AND. INTB.LE.153 .OR.
$ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 + 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 * ASCII is assumed, on Prime machines - ZCODE is the ASCII code
* plus 128 of either lower or upper case 'Z'. * plus 128 of either lower or upper case 'Z'.
* *
IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF (INTA.GE.225 .AND. INTA.LE.250) INTA = INTA - 32
IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 IF (INTB.GE.225 .AND. INTB.LE.250) INTB = INTB - 32
END IF END IF
LSAME = INTA.EQ.INTB LSAME = INTA .EQ. INTB
* *
* RETURN * RETURN
* *

97
lib/linalg/xerbla.f
View File

@ -1,34 +1,85 @@
*> \brief \b XERBLA
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download XERBLA + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/xerbla.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/xerbla.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/xerbla.f">
*> [TXT]</a>
*> \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 ) SUBROUTINE XERBLA( SRNAME, INFO )
* *
* -- LAPACK auxiliary routine (preliminary version) -- * -- LAPACK auxiliary routine (version 3.4.0) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * -- LAPACK is a software package provided by Univ. of Tennessee, --
* November 2006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER*(*) SRNAME CHARACTER*(*) SRNAME
INTEGER INFO 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 .. * .. Intrinsic Functions ..