1C***********************************************************************
2C LHS (Latin Hypercube Sampling) UNIX Library/Standalone.
3C Copyright (c) 2004, Sandia Corporation.  Under the terms of Contract
4C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
5C retains certain rights in this software.
6C
7C This software is distributed under the GNU Lesser General Public License.
8C For more information, see the README file in the LHS directory.
9C***********************************************************************
10C     Last change:  SLD  28 Mar 101    9:10 am
11C****************************************************************
12C SUBROUTINE POSDEF IS USED IN THE POSITIVE DEFINITE CHECK
13C OF THE CORRELATION MATRIX
14C
15!LHS_EXPORT_DEC ATTRIBUTES DLLEXPORT::POSDEF
16      SUBROUTINE POSDEF(ITEST)
17cc    ITEST returns number of iterations <= 20                          sld01
18cc    POSDEF is called from:  LHS                                       sld01
19C     Added LHS_ prefix to avoid clash with LAPACK - BMA-13
20cc    POSDEF is calls routines:  LHS_SSPEV,FINDIT                       sld01
21cc    POSDEF routine contains code for function PYTHAG                  sld01
22C     INCLUDE 'KILLFILE.INC'                                            GDW-96
23      USE KILLFILE
24C
25C     INCLUDE 'PARMS.INC'                                               GDW-96
26      USE PARMS
27cc    PARMS provides:  NVAR                                             sld01
28C     INCLUDE 'CCMATR.INC'                                              GDW-96
29      USE CCMATR
30cc    CCMATR provides:  CORR                                            sld01
31C
32C     These statements removed to make modules work - GDW-96
33C     COMMON/PDMAT/Z(NVAR,NVAR),D(NVAR)
34      USE PDMAT
35cc    PDMAT provides: Z and D arrays                                    sld01
36c
37      USE LOCALVARS, ONLY: WK
38cc    LOCALVARS provides:  WK array                                     sld01
39C
40      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
41C     DIMENSION WK((NVAR*(NVAR+1))/2)
42c     Moved to the LOCALVARS module
43C      DOUBLE PRECISION, ALLOCATABLE :: WK(:)
44C     NN = THE DIMENSIONS OF Z,  I.E.  Z(NN,NN).
45C     DATA NN/NVAR/
46C
47C
48C     EIG = THE VALUE THAT THE NEGATIVE EIGENVALUES ARE SET EQUAL TO.
49      DATA EIG/0.001/
50C
51C       These statements added to make modules work - GDW-96
52C       Allocate and initialize this work array
53c       Moved to the LOCALVARS module
54C        ALLOCATE( WK((NVAR*(NVAR+1))/2) )
55C        WK = 0.0
56C       NN = THE DIMENSIONS OF Z,  I.E.  Z(NN,NN).
57        NN = NVAR
58C
59C     M = MAXIMUM NUMBER OF ITERATIONS ALLOWED.
60      M=20
61C
62      NP=NCM
63c      NKX=(NP*(NP+1))/2       -----  removed 1-12-96, not used
64      ITEST=0
65      ICONV=0
66  100 CONTINUE
67      ITEST=ITEST+1
68      IF(ITEST.GT.M)THEN
69        WRITE(4,1000)
70        WRITE(99,1000)
71        KLLERR = .TRUE.
72        RETURN
73      ENDIF
74      REWIND 3
75      WRITE(3)CORR
76C     Added LHS_ prefix to avoid clash with LAPACK - BMA-13
77      CALL LHS_SSPEV(CORR,NP,D,Z,NN,WK,1,INFO)
78      If(KLLERR) Then
79         Return
80      END If
81      CALL FINDIT(NP,NN,EIG,ICONV)
82cc      If(KLLERR) Return -- FINDIT has no error conditions             sld01
83      IF(ICONV.EQ.0)GO TO 100
84      REWIND 3
85      READ(3)CORR
86C
87c      Moved to the LOCALVARS module
88C      DEALLOCATE( WK )
89C
90      RETURN
91 1000 FORMAT(1H1,'THE INPUT RANK CORRELATION MATRIX IS NOT POSITIVE ',
92     1       'DEFINITE.',/,' AN ITERATIVE PROCEDURE HAS FAILED TO ',
93     2       'PRODUCE A POSITIVE DEFINITE MATRIX AFTER 20 ITERATIONS.',
94     3       /,' THEREFORE, THE PROGRAM HAS BEEN TERMINATED.',/,' THE',
95     4       ' USER NEEDS TO REEVALUATE THE RELATIONSHIP OF THE ',
96     5       'CORRELATED VARIABLES IN THE MATRIX.')
97      END
98C****************************************************************
99C FUNCTION PYTHAG IS USED IN THE POSTIVE DEFINITE CHECK OF THE
100C CORRELATION MATRIX
101C
102      DOUBLE PRECISION FUNCTION PYTHAG(A,B)
103C***BEGIN PROLOGUE  PYTHAG
104C***REFER TO  EISDOC
105C     FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
106C***ROUTINES CALLED    (NONE)
107C***END PROLOGUE  PYTHAG
108      DOUBLE PRECISION A,B
109C
110      DOUBLE PRECISION P,Q,R,S,T
111C***FIRST EXECUTABLE STATEMENT  PYTHAG
112      P = MAX(ABS(A),ABS(B))
113      Q = MIN(ABS(A),ABS(B))
114      IF (Q .EQ. 0.0E0) GO TO 20
115   10 CONTINUE
116         R = (Q/P)**2
117         T = 4.0E0 + R
118         IF (T .EQ. 4.0E0) GO TO 20
119         S = R/T
120         P = P + 2.0E0*P*S
121         Q = Q*S
122      GO TO 10
123   20 PYTHAG = P
124      RETURN
125      END
126