← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/ddoglg.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/ddoglg.f
 
1
*DECK DDOGLG
 
2
      SUBROUTINE DDOGLG (N, R, LR, DIAG, QTB, DELTA, X, WA1, WA2)
 
3
C***BEGIN PROLOGUE  DDOGLG
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Subsidiary to DNSQ and DNSQE
 
6
C***LIBRARY   SLATEC
 
7
C***TYPE      DOUBLE PRECISION (DOGLEG-S, DDOGLG-D)
 
8
C***AUTHOR  (UNKNOWN)
 
9
C***DESCRIPTION
 
10
C
 
11
C     Given an M by N matrix A, an N by N nonsingular diagonal
 
12
C     matrix D, an M-vector B, and a positive number DELTA, the
 
13
C     problem is to determine the convex combination X of the
 
14
C     Gauss-Newton and scaled gradient directions that minimizes
 
15
C     (A*X - B) in the least squares sense, subject to the
 
16
C     restriction that the Euclidean norm of D*X be at most DELTA.
 
17
C
 
18
C     This subroutine completes the solution of the problem
 
19
C     if it is provided with the necessary information from the
 
20
C     QR factorization of A. That is, if A = Q*R, where Q has
 
21
C     orthogonal columns and R is an upper triangular matrix,
 
22
C     then DDOGLG expects the full upper triangle of R and
 
23
C     the first N components of (Q transpose)*B.
 
24
C
 
25
C     The subroutine statement is
 
26
C
 
27
C       SUBROUTINE DDOGLG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
 
28
C
 
29
C     where
 
30
C
 
31
C       N is a positive integer input variable set to the order of R.
 
32
C
 
33
C       R is an input array of length LR which must contain the upper
 
34
C         triangular matrix R stored by rows.
 
35
C
 
36
C       LR is a positive integer input variable not less than
 
37
C         (N*(N+1))/2.
 
38
C
 
39
C       DIAG is an input array of length N which must contain the
 
40
C         diagonal elements of the matrix D.
 
41
C
 
42
C       QTB is an input array of length N which must contain the first
 
43
C         N elements of the vector (Q transpose)*B.
 
44
C
 
45
C       DELTA is a positive input variable which specifies an upper
 
46
C         bound on the Euclidean norm of D*X.
 
47
C
 
48
C       X is an output array of length N which contains the desired
 
49
C         convex combination of the Gauss-Newton direction and the
 
50
C         scaled gradient direction.
 
51
C
 
52
C       WA1 and WA2 are work arrays of length N.
 
53
C
 
54
C***SEE ALSO  DNSQ, DNSQE
 
55
C***ROUTINES CALLED  D1MACH, DENORM
 
56
C***REVISION HISTORY  (YYMMDD)
 
57
C   800301  DATE WRITTEN
 
58
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
59
C   890831  Modified array declarations.  (WRB)
 
60
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
61
C   900326  Removed duplicate information from DESCRIPTION section.
 
62
C           (WRB)
 
63
C   900328  Added TYPE section.  (WRB)
 
64
C***END PROLOGUE  DDOGLG
 
65
      DOUBLE PRECISION D1MACH,DENORM
 
66
      INTEGER I, J, JJ, JP1, K, L, LR, N
 
67
      DOUBLE PRECISION ALPHA, BNORM, DELTA, DIAG(*), EPSMCH, GNORM,
 
68
     1     ONE, QNORM, QTB(*), R(*), SGNORM, SUM, TEMP, WA1(*),
 
69
     2     WA2(*), X(*), ZERO
 
70
      SAVE ONE, ZERO
 
71
      DATA ONE,ZERO /1.0D0,0.0D0/
 
72
C
 
73
C     EPSMCH IS THE MACHINE PRECISION.
 
74
C
 
75
C***FIRST EXECUTABLE STATEMENT  DDOGLG
 
76
      EPSMCH = D1MACH(4)
 
77
C
 
78
C     FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION.
 
79
C
 
80
      JJ = (N*(N + 1))/2 + 1
 
81
      DO 50 K = 1, N
 
82
         J = N - K + 1
 
83
         JP1 = J + 1
 
84
         JJ = JJ - K
 
85
         L = JJ + 1
 
86
         SUM = ZERO
 
87
         IF (N .LT. JP1) GO TO 20
 
88
         DO 10 I = JP1, N
 
89
            SUM = SUM + R(L)*X(I)
 
90
            L = L + 1
 
91
   10       CONTINUE
 
92
   20    CONTINUE
 
93
         TEMP = R(JJ)
 
94
         IF (TEMP .NE. ZERO) GO TO 40
 
95
         L = J
 
96
         DO 30 I = 1, J
 
97
            TEMP = MAX(TEMP,ABS(R(L)))
 
98
            L = L + N - I
 
99
   30       CONTINUE
 
100
         TEMP = EPSMCH*TEMP
 
101
         IF (TEMP .EQ. ZERO) TEMP = EPSMCH
 
102
   40    CONTINUE
 
103
         X(J) = (QTB(J) - SUM)/TEMP
 
104
   50    CONTINUE
 
105
C
 
106
C     TEST WHETHER THE GAUSS-NEWTON DIRECTION IS ACCEPTABLE.
 
107
C
 
108
      DO 60 J = 1, N
 
109
         WA1(J) = ZERO
 
110
         WA2(J) = DIAG(J)*X(J)
 
111
   60    CONTINUE
 
112
      QNORM = DENORM(N,WA2)
 
113
      IF (QNORM .LE. DELTA) GO TO 140
 
114
C
 
115
C     THE GAUSS-NEWTON DIRECTION IS NOT ACCEPTABLE.
 
116
C     NEXT, CALCULATE THE SCALED GRADIENT DIRECTION.
 
117
C
 
118
      L = 1
 
119
      DO 80 J = 1, N
 
120
         TEMP = QTB(J)
 
121
         DO 70 I = J, N
 
122
            WA1(I) = WA1(I) + R(L)*TEMP
 
123
            L = L + 1
 
124
   70       CONTINUE
 
125
         WA1(J) = WA1(J)/DIAG(J)
 
126
   80    CONTINUE
 
127
C
 
128
C     CALCULATE THE NORM OF THE SCALED GRADIENT AND TEST FOR
 
129
C     THE SPECIAL CASE IN WHICH THE SCALED GRADIENT IS ZERO.
 
130
C
 
131
      GNORM = DENORM(N,WA1)
 
132
      SGNORM = ZERO
 
133
      ALPHA = DELTA/QNORM
 
134
      IF (GNORM .EQ. ZERO) GO TO 120
 
135
C
 
136
C     CALCULATE THE POINT ALONG THE SCALED GRADIENT
 
137
C     AT WHICH THE QUADRATIC IS MINIMIZED.
 
138
C
 
139
      DO 90 J = 1, N
 
140
         WA1(J) = (WA1(J)/GNORM)/DIAG(J)
 
141
   90    CONTINUE
 
142
      L = 1
 
143
      DO 110 J = 1, N
 
144
         SUM = ZERO
 
145
         DO 100 I = J, N
 
146
            SUM = SUM + R(L)*WA1(I)
 
147
            L = L + 1
 
148
  100       CONTINUE
 
149
         WA2(J) = SUM
 
150
  110    CONTINUE
 
151
      TEMP = DENORM(N,WA2)
 
152
      SGNORM = (GNORM/TEMP)/TEMP
 
153
C
 
154
C     TEST WHETHER THE SCALED GRADIENT DIRECTION IS ACCEPTABLE.
 
155
C
 
156
      ALPHA = ZERO
 
157
      IF (SGNORM .GE. DELTA) GO TO 120
 
158
C
 
159
C     THE SCALED GRADIENT DIRECTION IS NOT ACCEPTABLE.
 
160
C     FINALLY, CALCULATE THE POINT ALONG THE DOGLEG
 
161
C     AT WHICH THE QUADRATIC IS MINIMIZED.
 
162
C
 
163
      BNORM = DENORM(N,QTB)
 
164
      TEMP = (BNORM/GNORM)*(BNORM/QNORM)*(SGNORM/DELTA)
 
165
      TEMP = TEMP - (DELTA/QNORM)*(SGNORM/DELTA)**2
 
166
     1       + SQRT((TEMP-(DELTA/QNORM))**2
 
167
     2               +(ONE-(DELTA/QNORM)**2)*(ONE-(SGNORM/DELTA)**2))
 
168
      ALPHA = ((DELTA/QNORM)*(ONE - (SGNORM/DELTA)**2))/TEMP
 
169
  120 CONTINUE
 
170
C
 
171
C     FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON
 
172
C     DIRECTION AND THE SCALED GRADIENT DIRECTION.
 
173
C
 
174
      TEMP = (ONE - ALPHA)*MIN(SGNORM,DELTA)
 
175
      DO 130 J = 1, N
 
176
         X(J) = TEMP*WA1(J) + ALPHA*X(J)
 
177
  130    CONTINUE
 
178
  140 CONTINUE
 
179
      RETURN
 
180
C
 
181
C     LAST CARD OF SUBROUTINE DDOGLG.
 
182
C
 
183
      END