← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dheqr.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/dheqr.f
 
1
*DECK DHEQR
 
2
      SUBROUTINE DHEQR (A, LDA, N, Q, INFO, IJOB)
 
3
C***BEGIN PROLOGUE  DHEQR
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Internal routine for DGMRES.
 
6
C***LIBRARY   SLATEC (SLAP)
 
7
C***CATEGORY  D2A4, D2B4
 
8
C***TYPE      DOUBLE PRECISION (SHEQR-S, DHEQR-D)
 
9
C***KEYWORDS  GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
 
10
C             NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
 
11
C***AUTHOR  Brown, Peter, (LLNL), pnbrown@llnl.gov
 
12
C           Hindmarsh, Alan, (LLNL), alanh@llnl.gov
 
13
C           Seager, Mark K., (LLNL), seager@llnl.gov
 
14
C             Lawrence Livermore National Laboratory
 
15
C             PO Box 808, L-60
 
16
C             Livermore, CA 94550 (510) 423-3141
 
17
C***DESCRIPTION
 
18
C        This   routine  performs  a QR   decomposition  of an  upper
 
19
C        Hessenberg matrix A using Givens  rotations.  There  are two
 
20
C        options  available: 1)  Performing  a fresh decomposition 2)
 
21
C        updating the QR factors by adding a row and  a column to the
 
22
C        matrix A.
 
23
C
 
24
C *Usage:
 
25
C      INTEGER LDA, N, INFO, IJOB
 
26
C      DOUBLE PRECISION A(LDA,N), Q(2*N)
 
27
C
 
28
C      CALL DHEQR(A, LDA, N, Q, INFO, IJOB)
 
29
C
 
30
C *Arguments:
 
31
C A      :INOUT    Double Precision A(LDA,N)
 
32
C         On input, the matrix to be decomposed.
 
33
C         On output, the upper triangular matrix R.
 
34
C         The factorization can be written Q*A = R, where
 
35
C         Q is a product of Givens rotations and R is upper
 
36
C         triangular.
 
37
C LDA    :IN       Integer
 
38
C         The leading dimension of the array A.
 
39
C N      :IN       Integer
 
40
C         A is an (N+1) by N Hessenberg matrix.
 
41
C Q      :OUT      Double Precision Q(2*N)
 
42
C         The factors c and s of each Givens rotation used
 
43
C         in decomposing A.
 
44
C INFO   :OUT      Integer
 
45
C         = 0  normal value.
 
46
C         = K  if  A(K,K) .eq. 0.0 .  This is not an error
 
47
C           condition for this subroutine, but it does
 
48
C           indicate that DHELS will divide by zero
 
49
C           if called.
 
50
C IJOB   :IN       Integer
 
51
C         = 1     means that a fresh decomposition of the
 
52
C                 matrix A is desired.
 
53
C         .ge. 2  means that the current decomposition of A
 
54
C                 will be updated by the addition of a row
 
55
C                 and a column.
 
56
C
 
57
C***SEE ALSO  DGMRES
 
58
C***ROUTINES CALLED  (NONE)
 
59
C***REVISION HISTORY  (YYMMDD)
 
60
C   890404  DATE WRITTEN
 
61
C   890404  Previous REVISION DATE
 
62
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
 
63
C   890922  Numerous changes to prologue to make closer to SLATEC
 
64
C           standard.  (FNF)
 
65
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
 
66
C   910411  Prologue converted to Version 4.0 format.  (BAB)
 
67
C   910506  Made subsidiary to DGMRES.  (FNF)
 
68
C   920511  Added complete declaration section.  (WRB)
 
69
C***END PROLOGUE  DHEQR
 
70
C         The following is for optimized compilation on LLNL/LTSS Crays.
 
71
CLLL. OPTIMIZE
 
72
C     .. Scalar Arguments ..
 
73
      INTEGER IJOB, INFO, LDA, N
 
74
C     .. Array Arguments ..
 
75
      DOUBLE PRECISION A(LDA,*), Q(*)
 
76
C     .. Local Scalars ..
 
77
      DOUBLE PRECISION C, S, T, T1, T2
 
78
      INTEGER I, IQ, J, K, KM1, KP1, NM1
 
79
C     .. Intrinsic Functions ..
 
80
      INTRINSIC ABS, SQRT
 
81
C***FIRST EXECUTABLE STATEMENT  DHEQR
 
82
      IF (IJOB .GT. 1) GO TO 70
 
83
C   -------------------------------------------------------------------
 
84
C         A new factorization is desired.
 
85
C   -------------------------------------------------------------------
 
86
C         QR decomposition without pivoting.
 
87
C
 
88
      INFO = 0
 
89
      DO 60 K = 1, N
 
90
         KM1 = K - 1
 
91
         KP1 = K + 1
 
92
C
 
93
C           Compute K-th column of R.
 
94
C           First, multiply the K-th column of A by the previous
 
95
C           K-1 Givens rotations.
 
96
C
 
97
         IF (KM1 .LT. 1) GO TO 20
 
98
         DO 10 J = 1, KM1
 
99
            I = 2*(J-1) + 1
 
100
            T1 = A(J,K)
 
101
            T2 = A(J+1,K)
 
102
            C = Q(I)
 
103
            S = Q(I+1)
 
104
            A(J,K) = C*T1 - S*T2
 
105
            A(J+1,K) = S*T1 + C*T2
 
106
 10      CONTINUE
 
107
C
 
108
C         Compute Givens components C and S.
 
109
C
 
110
 20      CONTINUE
 
111
         IQ = 2*KM1 + 1
 
112
         T1 = A(K,K)
 
113
         T2 = A(KP1,K)
 
114
         IF( T2.EQ.0.0D0 ) THEN
 
115
            C = 1
 
116
            S = 0
 
117
         ELSEIF( ABS(T2).GE.ABS(T1) ) THEN
 
118
            T = T1/T2
 
119
            S = -1.0D0/SQRT(1.0D0+T*T)
 
120
            C = -S*T
 
121
         ELSE
 
122
            T = T2/T1
 
123
            C = 1.0D0/SQRT(1.0D0+T*T)
 
124
            S = -C*T
 
125
         ENDIF
 
126
         Q(IQ) = C
 
127
         Q(IQ+1) = S
 
128
         A(K,K) = C*T1 - S*T2
 
129
         IF( A(K,K).EQ.0.0D0 ) INFO = K
 
130
 60   CONTINUE
 
131
      RETURN
 
132
C   -------------------------------------------------------------------
 
133
C         The old factorization of a will be updated.  A row and a
 
134
C         column has been added to the matrix A.  N by N-1 is now
 
135
C         the old size of the matrix.
 
136
C   -------------------------------------------------------------------
 
137
 70   CONTINUE
 
138
      NM1 = N - 1
 
139
C   -------------------------------------------------------------------
 
140
C         Multiply the new column by the N previous Givens rotations.
 
141
C   -------------------------------------------------------------------
 
142
      DO 100 K = 1,NM1
 
143
         I = 2*(K-1) + 1
 
144
         T1 = A(K,N)
 
145
         T2 = A(K+1,N)
 
146
         C = Q(I)
 
147
         S = Q(I+1)
 
148
         A(K,N) = C*T1 - S*T2
 
149
         A(K+1,N) = S*T1 + C*T2
 
150
 100  CONTINUE
 
151
C   -------------------------------------------------------------------
 
152
C         Complete update of decomposition by forming last Givens
 
153
C         rotation, and multiplying it times the column
 
154
C         vector(A(N,N),A(NP1,N)).
 
155
C   -------------------------------------------------------------------
 
156
      INFO = 0
 
157
      T1 = A(N,N)
 
158
      T2 = A(N+1,N)
 
159
      IF ( T2.EQ.0.0D0 ) THEN
 
160
         C = 1
 
161
         S = 0
 
162
      ELSEIF( ABS(T2).GE.ABS(T1) ) THEN
 
163
         T = T1/T2
 
164
         S = -1.0D0/SQRT(1.0D0+T*T)
 
165
         C = -S*T
 
166
      ELSE
 
167
         T = T2/T1
 
168
         C = 1.0D0/SQRT(1.0D0+T*T)
 
169
         S = -C*T
 
170
      ENDIF
 
171
      IQ = 2*N - 1
 
172
      Q(IQ) = C
 
173
      Q(IQ+1) = S
 
174
      A(N,N) = C*T1 - S*T2
 
175
      IF (A(N,N) .EQ. 0.0D0) INFO = N
 
176
      RETURN
 
177
C------------- LAST LINE OF DHEQR FOLLOWS ----------------------------
 
178
      END