← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dsisl.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/dsisl.f
 
1
*DECK DSISL
 
2
      SUBROUTINE DSISL (A, LDA, N, KPVT, B)
 
3
C***BEGIN PROLOGUE  DSISL
 
4
C***PURPOSE  Solve a real symmetric system using the factors obtained
 
5
C            from SSIFA.
 
6
C***LIBRARY   SLATEC (LINPACK)
 
7
C***CATEGORY  D2B1A
 
8
C***TYPE      DOUBLE PRECISION (SSISL-S, DSISL-D, CHISL-C, CSISL-C)
 
9
C***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE, SYMMETRIC
 
10
C***AUTHOR  Bunch, J., (UCSD)
 
11
C***DESCRIPTION
 
12
C
 
13
C     DSISL solves the double precision symmetric system
 
14
C     A * X = B
 
15
C     using the factors computed by DSIFA.
 
16
C
 
17
C     On Entry
 
18
C
 
19
C        A       DOUBLE PRECISION(LDA,N)
 
20
C                the output from DSIFA.
 
21
C
 
22
C        LDA     INTEGER
 
23
C                the leading dimension of the array  A .
 
24
C
 
25
C        N       INTEGER
 
26
C                the order of the matrix  A .
 
27
C
 
28
C        KPVT    INTEGER(N)
 
29
C                the pivot vector from DSIFA.
 
30
C
 
31
C        B       DOUBLE PRECISION(N)
 
32
C                the right hand side vector.
 
33
C
 
34
C     On Return
 
35
C
 
36
C        B       the solution vector  X .
 
37
C
 
38
C     Error Condition
 
39
C
 
40
C        A division by zero may occur if  DSICO  has set RCOND .EQ. 0.0
 
41
C        or  DSIFA  has set INFO .NE. 0  .
 
42
C
 
43
C     To compute  INVERSE(A) * C  where  C  is a matrix
 
44
C     with  P  columns
 
45
C           CALL DSIFA(A,LDA,N,KPVT,INFO)
 
46
C           IF (INFO .NE. 0) GO TO ...
 
47
C           DO 10 J = 1, P
 
48
C              CALL DSISL(A,LDA,N,KPVT,C(1,J))
 
49
C        10 CONTINUE
 
50
C
 
51
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
52
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
53
C***ROUTINES CALLED  DAXPY, DDOT
 
54
C***REVISION HISTORY  (YYMMDD)
 
55
C   780814  DATE WRITTEN
 
56
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
57
C   890831  Modified array declarations.  (WRB)
 
58
C   891107  Modified routine equivalence list.  (WRB)
 
59
C   891107  REVISION DATE from Version 3.2
 
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   920501  Reformatted the REFERENCES section.  (WRB)
 
64
C***END PROLOGUE  DSISL
 
65
      INTEGER LDA,N,KPVT(*)
 
66
      DOUBLE PRECISION A(LDA,*),B(*)
 
67
C
 
68
      DOUBLE PRECISION AK,AKM1,BK,BKM1,DDOT,DENOM,TEMP
 
69
      INTEGER K,KP
 
70
C
 
71
C     LOOP BACKWARD APPLYING THE TRANSFORMATIONS AND
 
72
C     D INVERSE TO B.
 
73
C
 
74
C***FIRST EXECUTABLE STATEMENT  DSISL
 
75
      K = N
 
76
   10 IF (K .EQ. 0) GO TO 80
 
77
         IF (KPVT(K) .LT. 0) GO TO 40
 
78
C
 
79
C           1 X 1 PIVOT BLOCK.
 
80
C
 
81
            IF (K .EQ. 1) GO TO 30
 
82
               KP = KPVT(K)
 
83
               IF (KP .EQ. K) GO TO 20
 
84
C
 
85
C                 INTERCHANGE.
 
86
C
 
87
                  TEMP = B(K)
 
88
                  B(K) = B(KP)
 
89
                  B(KP) = TEMP
 
90
   20          CONTINUE
 
91
C
 
92
C              APPLY THE TRANSFORMATION.
 
93
C
 
94
               CALL DAXPY(K-1,B(K),A(1,K),1,B(1),1)
 
95
   30       CONTINUE
 
96
C
 
97
C           APPLY D INVERSE.
 
98
C
 
99
            B(K) = B(K)/A(K,K)
 
100
            K = K - 1
 
101
         GO TO 70
 
102
   40    CONTINUE
 
103
C
 
104
C           2 X 2 PIVOT BLOCK.
 
105
C
 
106
            IF (K .EQ. 2) GO TO 60
 
107
               KP = ABS(KPVT(K))
 
108
               IF (KP .EQ. K - 1) GO TO 50
 
109
C
 
110
C                 INTERCHANGE.
 
111
C
 
112
                  TEMP = B(K-1)
 
113
                  B(K-1) = B(KP)
 
114
                  B(KP) = TEMP
 
115
   50          CONTINUE
 
116
C
 
117
C              APPLY THE TRANSFORMATION.
 
118
C
 
119
               CALL DAXPY(K-2,B(K),A(1,K),1,B(1),1)
 
120
               CALL DAXPY(K-2,B(K-1),A(1,K-1),1,B(1),1)
 
121
   60       CONTINUE
 
122
C
 
123
C           APPLY D INVERSE.
 
124
C
 
125
            AK = A(K,K)/A(K-1,K)
 
126
            AKM1 = A(K-1,K-1)/A(K-1,K)
 
127
            BK = B(K)/A(K-1,K)
 
128
            BKM1 = B(K-1)/A(K-1,K)
 
129
            DENOM = AK*AKM1 - 1.0D0
 
130
            B(K) = (AKM1*BK - BKM1)/DENOM
 
131
            B(K-1) = (AK*BKM1 - BK)/DENOM
 
132
            K = K - 2
 
133
   70    CONTINUE
 
134
      GO TO 10
 
135
   80 CONTINUE
 
136
C
 
137
C     LOOP FORWARD APPLYING THE TRANSFORMATIONS.
 
138
C
 
139
      K = 1
 
140
   90 IF (K .GT. N) GO TO 160
 
141
         IF (KPVT(K) .LT. 0) GO TO 120
 
142
C
 
143
C           1 X 1 PIVOT BLOCK.
 
144
C
 
145
            IF (K .EQ. 1) GO TO 110
 
146
C
 
147
C              APPLY THE TRANSFORMATION.
 
148
C
 
149
               B(K) = B(K) + DDOT(K-1,A(1,K),1,B(1),1)
 
150
               KP = KPVT(K)
 
151
               IF (KP .EQ. K) GO TO 100
 
152
C
 
153
C                 INTERCHANGE.
 
154
C
 
155
                  TEMP = B(K)
 
156
                  B(K) = B(KP)
 
157
                  B(KP) = TEMP
 
158
  100          CONTINUE
 
159
  110       CONTINUE
 
160
            K = K + 1
 
161
         GO TO 150
 
162
  120    CONTINUE
 
163
C
 
164
C           2 X 2 PIVOT BLOCK.
 
165
C
 
166
            IF (K .EQ. 1) GO TO 140
 
167
C
 
168
C              APPLY THE TRANSFORMATION.
 
169
C
 
170
               B(K) = B(K) + DDOT(K-1,A(1,K),1,B(1),1)
 
171
               B(K+1) = B(K+1) + DDOT(K-1,A(1,K+1),1,B(1),1)
 
172
               KP = ABS(KPVT(K))
 
173
               IF (KP .EQ. K) GO TO 130
 
174
C
 
175
C                 INTERCHANGE.
 
176
C
 
177
                  TEMP = B(K)
 
178
                  B(K) = B(KP)
 
179
                  B(KP) = TEMP
 
180
  130          CONTINUE
 
181
  140       CONTINUE
 
182
            K = K + 2
 
183
  150    CONTINUE
 
184
      GO TO 90
 
185
  160 CONTINUE
 
186
      RETURN
 
187
      END