added lapack 3.7.0 with latest patches from git
This commit is contained in:
Werner Saar
@@ -0,0 +1,278 @@
|
||||
*> \brief \b CLAHILB
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
*
|
||||
* Online html documentation available at
|
||||
* http://www.netlib.org/lapack/explore-html/
|
||||
*
|
||||
* Definition:
|
||||
* ===========
|
||||
*
|
||||
* SUBROUTINE CLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
|
||||
* INFO, PATH)
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
* INTEGER N, NRHS, LDA, LDX, LDB, INFO
|
||||
* .. Array Arguments ..
|
||||
* REAL WORK(N)
|
||||
* COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
|
||||
* CHARACTER*3 PATH
|
||||
* ..
|
||||
*
|
||||
*
|
||||
*> \par Purpose:
|
||||
* =============
|
||||
*>
|
||||
*> \verbatim
|
||||
*>
|
||||
*> CLAHILB generates an N by N scaled Hilbert matrix in A along with
|
||||
*> NRHS right-hand sides in B and solutions in X such that A*X=B.
|
||||
*>
|
||||
*> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
|
||||
*> entries are integers. The right-hand sides are the first NRHS
|
||||
*> columns of M * the identity matrix, and the solutions are the
|
||||
*> first NRHS columns of the inverse Hilbert matrix.
|
||||
*>
|
||||
*> The condition number of the Hilbert matrix grows exponentially with
|
||||
*> its size, roughly as O(e ** (3.5*N)). Additionally, the inverse
|
||||
*> Hilbert matrices beyond a relatively small dimension cannot be
|
||||
*> generated exactly without extra precision. Precision is exhausted
|
||||
*> when the largest entry in the inverse Hilbert matrix is greater than
|
||||
*> 2 to the power of the number of bits in the fraction of the data type
|
||||
*> used plus one, which is 24 for single precision.
|
||||
*>
|
||||
*> In single, the generated solution is exact for N <= 6 and has
|
||||
*> small componentwise error for 7 <= N <= 11.
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
* ==========
|
||||
*
|
||||
*> \param[in] N
|
||||
*> \verbatim
|
||||
*> N is INTEGER
|
||||
*> The dimension of the matrix A.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] NRHS
|
||||
*> \verbatim
|
||||
*> NRHS is INTEGER
|
||||
*> The requested number of right-hand sides.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[out] A
|
||||
*> \verbatim
|
||||
*> A is COMPLEX array, dimension (LDA, N)
|
||||
*> The generated scaled Hilbert matrix.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] LDA
|
||||
*> \verbatim
|
||||
*> LDA is INTEGER
|
||||
*> The leading dimension of the array A. LDA >= N.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[out] X
|
||||
*> \verbatim
|
||||
*> X is COMPLEX array, dimension (LDX, NRHS)
|
||||
*> The generated exact solutions. Currently, the first NRHS
|
||||
*> columns of the inverse Hilbert matrix.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] LDX
|
||||
*> \verbatim
|
||||
*> LDX is INTEGER
|
||||
*> The leading dimension of the array X. LDX >= N.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[out] B
|
||||
*> \verbatim
|
||||
*> B is REAL array, dimension (LDB, NRHS)
|
||||
*> The generated right-hand sides. Currently, the first NRHS
|
||||
*> columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] LDB
|
||||
*> \verbatim
|
||||
*> LDB is INTEGER
|
||||
*> The leading dimension of the array B. LDB >= N.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[out] WORK
|
||||
*> \verbatim
|
||||
*> WORK is REAL array, dimension (N)
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[out] INFO
|
||||
*> \verbatim
|
||||
*> INFO is INTEGER
|
||||
*> = 0: successful exit
|
||||
*> = 1: N is too large; the data is still generated but may not
|
||||
*> be not exact.
|
||||
*> < 0: if INFO = -i, the i-th argument had an illegal value
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] PATH
|
||||
*> \verbatim
|
||||
*> PATH is CHARACTER*3
|
||||
*> The LAPACK path name.
|
||||
*> \endverbatim
|
||||
*
|
||||
* Authors:
|
||||
* ========
|
||||
*
|
||||
*> \author Univ. of Tennessee
|
||||
*> \author Univ. of California Berkeley
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date December 2016
|
||||
*
|
||||
*> \ingroup complex_matgen
|
||||
*
|
||||
* =====================================================================
|
||||
SUBROUTINE CLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
|
||||
$ INFO, PATH)
|
||||
*
|
||||
* -- LAPACK test routine (version 3.7.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* December 2016
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER N, NRHS, LDA, LDX, LDB, INFO
|
||||
* .. Array Arguments ..
|
||||
REAL WORK(N)
|
||||
COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
|
||||
CHARACTER*3 PATH
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
* .. Local Scalars ..
|
||||
INTEGER TM, TI, R
|
||||
INTEGER M
|
||||
INTEGER I, J
|
||||
COMPLEX TMP
|
||||
CHARACTER*2 C2
|
||||
|
||||
* .. Parameters ..
|
||||
* NMAX_EXACT the largest dimension where the generated data is
|
||||
* exact.
|
||||
* NMAX_APPROX the largest dimension where the generated data has
|
||||
* a small componentwise relative error.
|
||||
* ??? complex uses how many bits ???
|
||||
INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
|
||||
PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8)
|
||||
|
||||
* d's are generated from random permuation of those eight elements.
|
||||
COMPLEX D1(8), D2(8), INVD1(8), INVD2(8)
|
||||
DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
|
||||
DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
|
||||
|
||||
DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
|
||||
$ (-.5,-.5),(.5,-.5),(.5,.5)/
|
||||
DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
|
||||
$ (-.5,.5),(.5,.5),(.5,-.5)/
|
||||
|
||||
* ..
|
||||
* .. External Functions
|
||||
EXTERNAL CLASET, LSAMEN
|
||||
INTRINSIC REAL
|
||||
LOGICAL LSAMEN
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
C2 = PATH( 2: 3 )
|
||||
*
|
||||
* Test the input arguments
|
||||
*
|
||||
INFO = 0
|
||||
IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
|
||||
INFO = -1
|
||||
ELSE IF (NRHS .LT. 0) THEN
|
||||
INFO = -2
|
||||
ELSE IF (LDA .LT. N) THEN
|
||||
INFO = -4
|
||||
ELSE IF (LDX .LT. N) THEN
|
||||
INFO = -6
|
||||
ELSE IF (LDB .LT. N) THEN
|
||||
INFO = -8
|
||||
END IF
|
||||
IF (INFO .LT. 0) THEN
|
||||
CALL XERBLA('CLAHILB', -INFO)
|
||||
RETURN
|
||||
END IF
|
||||
IF (N .GT. NMAX_EXACT) THEN
|
||||
INFO = 1
|
||||
END IF
|
||||
|
||||
* Compute M = the LCM of the integers [1, 2*N-1]. The largest
|
||||
* reasonable N is small enough that integers suffice (up to N = 11).
|
||||
M = 1
|
||||
DO I = 2, (2*N-1)
|
||||
TM = M
|
||||
TI = I
|
||||
R = MOD(TM, TI)
|
||||
DO WHILE (R .NE. 0)
|
||||
TM = TI
|
||||
TI = R
|
||||
R = MOD(TM, TI)
|
||||
END DO
|
||||
M = (M / TI) * I
|
||||
END DO
|
||||
|
||||
* Generate the scaled Hilbert matrix in A
|
||||
* If we are testing SY routines, take
|
||||
* D1_i = D2_i, else, D1_i = D2_i*
|
||||
IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
|
||||
DO J = 1, N
|
||||
DO I = 1, N
|
||||
A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
|
||||
$ * D1(MOD(I,SIZE_D)+1)
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
DO J = 1, N
|
||||
DO I = 1, N
|
||||
A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
|
||||
$ * D2(MOD(I,SIZE_D)+1)
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
|
||||
* Generate matrix B as simply the first NRHS columns of M * the
|
||||
* identity.
|
||||
TMP = REAL(M)
|
||||
CALL CLASET('Full', N, NRHS, (0.0,0.0), TMP, B, LDB)
|
||||
|
||||
* Generate the true solutions in X. Because B = the first NRHS
|
||||
* columns of M*I, the true solutions are just the first NRHS columns
|
||||
* of the inverse Hilbert matrix.
|
||||
WORK(1) = N
|
||||
DO J = 2, N
|
||||
WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) )
|
||||
$ * (N +J -1)
|
||||
END DO
|
||||
|
||||
* If we are testing SY routines,
|
||||
* take D1_i = D2_i, else, D1_i = D2_i*
|
||||
IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
|
||||
DO J = 1, NRHS
|
||||
DO I = 1, N
|
||||
X(I, J) =
|
||||
$ INVD1(MOD(J,SIZE_D)+1) *
|
||||
$ ((WORK(I)*WORK(J)) / (I + J - 1))
|
||||
$ * INVD1(MOD(I,SIZE_D)+1)
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
DO J = 1, NRHS
|
||||
DO I = 1, N
|
||||
X(I, J) =
|
||||
$ INVD2(MOD(J,SIZE_D)+1) *
|
||||
$ ((WORK(I)*WORK(J)) / (I + J - 1))
|
||||
$ * INVD1(MOD(I,SIZE_D)+1)
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
END
|
||||
|
||||
Reference in New Issue
Block a user