← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dlpdp.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/dlpdp.f
 
1
*DECK DLPDP
 
2
      SUBROUTINE DLPDP (A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS,
 
3
     +   IS)
 
4
C***BEGIN PROLOGUE  DLPDP
 
5
C***SUBSIDIARY
 
6
C***PURPOSE  Subsidiary to DLSEI
 
7
C***LIBRARY   SLATEC
 
8
C***TYPE      DOUBLE PRECISION (LPDP-S, DLPDP-D)
 
9
C***AUTHOR  Hanson, R. J., (SNLA)
 
10
C           Haskell, K. H., (SNLA)
 
11
C***DESCRIPTION
 
12
C
 
13
C  **** Double Precision version of LPDP ****
 
14
C     DIMENSION A(MDA,N+1),PRGOPT(*),X(N),WS((M+2)*(N+7)),IS(M+N+1),
 
15
C     where N=N1+N2.  This is a slight overestimate for WS(*).
 
16
C
 
17
C     Determine an N1-vector W, and
 
18
C               an N2-vector Z
 
19
C     which minimizes the Euclidean length of W
 
20
C     subject to G*W+H*Z .GE. Y.
 
21
C     This is the least projected distance problem, LPDP.
 
22
C     The matrices G and H are of respective
 
23
C     dimensions M by N1 and M by N2.
 
24
C
 
25
C     Called by subprogram DLSI( ).
 
26
C
 
27
C     The matrix
 
28
C                (G H Y)
 
29
C
 
30
C     occupies rows 1,...,M and cols 1,...,N1+N2+1 of A(*,*).
 
31
C
 
32
C     The solution (W) is returned in X(*).
 
33
C                  (Z)
 
34
C
 
35
C     The value of MODE indicates the status of
 
36
C     the computation after returning to the user.
 
37
C
 
38
C          MODE=1  The solution was successfully obtained.
 
39
C
 
40
C          MODE=2  The inequalities are inconsistent.
 
41
C
 
42
C***SEE ALSO  DLSEI
 
43
C***ROUTINES CALLED  DCOPY, DDOT, DNRM2, DSCAL, DWNNLS
 
44
C***REVISION HISTORY  (YYMMDD)
 
45
C   790701  DATE WRITTEN
 
46
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
47
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
48
C   900328  Added TYPE section.  (WRB)
 
49
C   910408  Updated the AUTHOR section.  (WRB)
 
50
C***END PROLOGUE  DLPDP
 
51
C
 
52
      INTEGER I, IS(*), IW, IX, J, L, M, MDA, MODE, MODEW, N, N1, N2,
 
53
     *     NP1
 
54
      DOUBLE PRECISION A(MDA,*), DDOT, DNRM2, FAC, ONE,
 
55
     *     PRGOPT(*), RNORM, SC, WNORM, WS(*), X(*), YNORM, ZERO
 
56
      SAVE ZERO, ONE, FAC
 
57
      DATA ZERO,ONE /0.0D0,1.0D0/, FAC /0.1D0/
 
58
C***FIRST EXECUTABLE STATEMENT  DLPDP
 
59
      N = N1 + N2
 
60
      MODE = 1
 
61
      IF (M .GT. 0) GO TO 20
 
62
         IF (N .LE. 0) GO TO 10
 
63
            X(1) = ZERO
 
64
            CALL DCOPY(N,X,0,X,1)
 
65
   10    CONTINUE
 
66
         WNORM = ZERO
 
67
      GO TO 200
 
68
   20 CONTINUE
 
69
C        BEGIN BLOCK PERMITTING ...EXITS TO 190
 
70
            NP1 = N + 1
 
71
C
 
72
C           SCALE NONZERO ROWS OF INEQUALITY MATRIX TO HAVE LENGTH ONE.
 
73
            DO 40 I = 1, M
 
74
               SC = DNRM2(N,A(I,1),MDA)
 
75
               IF (SC .EQ. ZERO) GO TO 30
 
76
                  SC = ONE/SC
 
77
                  CALL DSCAL(NP1,SC,A(I,1),MDA)
 
78
   30          CONTINUE
 
79
   40       CONTINUE
 
80
C
 
81
C           SCALE RT.-SIDE VECTOR TO HAVE LENGTH ONE (OR ZERO).
 
82
            YNORM = DNRM2(M,A(1,NP1),1)
 
83
            IF (YNORM .EQ. ZERO) GO TO 50
 
84
               SC = ONE/YNORM
 
85
               CALL DSCAL(M,SC,A(1,NP1),1)
 
86
   50       CONTINUE
 
87
C
 
88
C           SCALE COLS OF MATRIX H.
 
89
            J = N1 + 1
 
90
   60       IF (J .GT. N) GO TO 70
 
91
               SC = DNRM2(M,A(1,J),1)
 
92
               IF (SC .NE. ZERO) SC = ONE/SC
 
93
               CALL DSCAL(M,SC,A(1,J),1)
 
94
               X(J) = SC
 
95
               J = J + 1
 
96
            GO TO 60
 
97
   70       CONTINUE
 
98
            IF (N1 .LE. 0) GO TO 130
 
99
C
 
100
C              COPY TRANSPOSE OF (H G Y) TO WORK ARRAY WS(*).
 
101
               IW = 0
 
102
               DO 80 I = 1, M
 
103
C
 
104
C                 MOVE COL OF TRANSPOSE OF H INTO WORK ARRAY.
 
105
                  CALL DCOPY(N2,A(I,N1+1),MDA,WS(IW+1),1)
 
106
                  IW = IW + N2
 
107
C
 
108
C                 MOVE COL OF TRANSPOSE OF G INTO WORK ARRAY.
 
109
                  CALL DCOPY(N1,A(I,1),MDA,WS(IW+1),1)
 
110
                  IW = IW + N1
 
111
C
 
112
C                 MOVE COMPONENT OF VECTOR Y INTO WORK ARRAY.
 
113
                  WS(IW+1) = A(I,NP1)
 
114
                  IW = IW + 1
 
115
   80          CONTINUE
 
116
               WS(IW+1) = ZERO
 
117
               CALL DCOPY(N,WS(IW+1),0,WS(IW+1),1)
 
118
               IW = IW + N
 
119
               WS(IW+1) = ONE
 
120
               IW = IW + 1
 
121
C
 
122
C              SOLVE EU=F SUBJECT TO (TRANSPOSE OF H)U=0, U.GE.0.  THE
 
123
C              MATRIX E = TRANSPOSE OF (G Y), AND THE (N+1)-VECTOR
 
124
C              F = TRANSPOSE OF (0,...,0,1).
 
125
               IX = IW + 1
 
126
               IW = IW + M
 
127
C
 
128
C              DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF
 
129
C              DWNNLS( ).
 
130
               IS(1) = 0
 
131
               IS(2) = 0
 
132
               CALL DWNNLS(WS,NP1,N2,NP1-N2,M,0,PRGOPT,WS(IX),RNORM,
 
133
     *                     MODEW,IS,WS(IW+1))
 
134
C
 
135
C              COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY W.
 
136
               SC = ONE - DDOT(M,A(1,NP1),1,WS(IX),1)
 
137
               IF (ONE + FAC*ABS(SC) .EQ. ONE .OR. RNORM .LE. ZERO)
 
138
     *            GO TO 110
 
139
                  SC = ONE/SC
 
140
                  DO 90 J = 1, N1
 
141
                     X(J) = SC*DDOT(M,A(1,J),1,WS(IX),1)
 
142
   90             CONTINUE
 
143
C
 
144
C                 COMPUTE THE VECTOR Q=Y-GW.  OVERWRITE Y WITH THIS
 
145
C                 VECTOR.
 
146
                  DO 100 I = 1, M
 
147
                     A(I,NP1) = A(I,NP1) - DDOT(N1,A(I,1),MDA,X,1)
 
148
  100             CONTINUE
 
149
               GO TO 120
 
150
  110          CONTINUE
 
151
                  MODE = 2
 
152
C        .........EXIT
 
153
                  GO TO 190
 
154
  120          CONTINUE
 
155
  130       CONTINUE
 
156
            IF (N2 .LE. 0) GO TO 180
 
157
C
 
158
C              COPY TRANSPOSE OF (H Q) TO WORK ARRAY WS(*).
 
159
               IW = 0
 
160
               DO 140 I = 1, M
 
161
                  CALL DCOPY(N2,A(I,N1+1),MDA,WS(IW+1),1)
 
162
                  IW = IW + N2
 
163
                  WS(IW+1) = A(I,NP1)
 
164
                  IW = IW + 1
 
165
  140          CONTINUE
 
166
               WS(IW+1) = ZERO
 
167
               CALL DCOPY(N2,WS(IW+1),0,WS(IW+1),1)
 
168
               IW = IW + N2
 
169
               WS(IW+1) = ONE
 
170
               IW = IW + 1
 
171
               IX = IW + 1
 
172
               IW = IW + M
 
173
C
 
174
C              SOLVE RV=S SUBJECT TO V.GE.0.  THE MATRIX R =(TRANSPOSE
 
175
C              OF (H Q)), WHERE Q=Y-GW.  THE (N2+1)-VECTOR S =(TRANSPOSE
 
176
C              OF (0,...,0,1)).
 
177
C
 
178
C              DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF
 
179
C              DWNNLS( ).
 
180
               IS(1) = 0
 
181
               IS(2) = 0
 
182
               CALL DWNNLS(WS,N2+1,0,N2+1,M,0,PRGOPT,WS(IX),RNORM,MODEW,
 
183
     *                     IS,WS(IW+1))
 
184
C
 
185
C              COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY Z.
 
186
               SC = ONE - DDOT(M,A(1,NP1),1,WS(IX),1)
 
187
               IF (ONE + FAC*ABS(SC) .EQ. ONE .OR. RNORM .LE. ZERO)
 
188
     *            GO TO 160
 
189
                  SC = ONE/SC
 
190
                  DO 150 J = 1, N2
 
191
                     L = N1 + J
 
192
                     X(L) = SC*DDOT(M,A(1,L),1,WS(IX),1)*X(L)
 
193
  150             CONTINUE
 
194
               GO TO 170
 
195
  160          CONTINUE
 
196
                  MODE = 2
 
197
C        .........EXIT
 
198
                  GO TO 190
 
199
  170          CONTINUE
 
200
  180       CONTINUE
 
201
C
 
202
C           ACCOUNT FOR SCALING OF RT.-SIDE VECTOR IN SOLUTION.
 
203
            CALL DSCAL(N,YNORM,X,1)
 
204
            WNORM = DNRM2(N1,X,1)
 
205
  190    CONTINUE
 
206
  200 CONTINUE
 
207
      RETURN
 
208
      END