← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dtrsl.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/dtrsl.f
 
1
*DECK DTRSL
 
2
      SUBROUTINE DTRSL (T, LDT, N, B, JOB, INFO)
 
3
C***BEGIN PROLOGUE  DTRSL
 
4
C***PURPOSE  Solve a system of the form  T*X=B or TRANS(T)*X=B, where
 
5
C            T is a triangular matrix.
 
6
C***LIBRARY   SLATEC (LINPACK)
 
7
C***CATEGORY  D2A3
 
8
C***TYPE      DOUBLE PRECISION (STRSL-S, DTRSL-D, CTRSL-C)
 
9
C***KEYWORDS  LINEAR ALGEBRA, LINPACK, TRIANGULAR LINEAR SYSTEM,
 
10
C             TRIANGULAR MATRIX
 
11
C***AUTHOR  Stewart, G. W., (U. of Maryland)
 
12
C***DESCRIPTION
 
13
C
 
14
C     DTRSL solves systems of the form
 
15
C
 
16
C                   T * X = B
 
17
C     or
 
18
C                   TRANS(T) * X = B
 
19
C
 
20
C     where T is a triangular matrix of order N.  Here TRANS(T)
 
21
C     denotes the transpose of the matrix T.
 
22
C
 
23
C     On Entry
 
24
C
 
25
C         T         DOUBLE PRECISION(LDT,N)
 
26
C                   T contains the matrix of the system.  The zero
 
27
C                   elements of the matrix are not referenced, and
 
28
C                   the corresponding elements of the array can be
 
29
C                   used to store other information.
 
30
C
 
31
C         LDT       INTEGER
 
32
C                   LDT is the leading dimension of the array T.
 
33
C
 
34
C         N         INTEGER
 
35
C                   N is the order of the system.
 
36
C
 
37
C         B         DOUBLE PRECISION(N).
 
38
C                   B contains the right hand side of the system.
 
39
C
 
40
C         JOB       INTEGER
 
41
C                   JOB specifies what kind of system is to be solved.
 
42
C                   If JOB is
 
43
C
 
44
C                        00   solve T*X=B, T lower triangular,
 
45
C                        01   solve T*X=B, T upper triangular,
 
46
C                        10   solve TRANS(T)*X=B, T lower triangular,
 
47
C                        11   solve TRANS(T)*X=B, T upper triangular.
 
48
C
 
49
C     On Return
 
50
C
 
51
C         B         B contains the solution, if INFO .EQ. 0.
 
52
C                   Otherwise B is unaltered.
 
53
C
 
54
C         INFO      INTEGER
 
55
C                   INFO contains zero if the system is nonsingular.
 
56
C                   Otherwise INFO contains the index of
 
57
C                   the first zero diagonal element of T.
 
58
C
 
59
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
60
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
61
C***ROUTINES CALLED  DAXPY, DDOT
 
62
C***REVISION HISTORY  (YYMMDD)
 
63
C   780814  DATE WRITTEN
 
64
C   890831  Modified array declarations.  (WRB)
 
65
C   890831  REVISION DATE from Version 3.2
 
66
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
67
C   900326  Removed duplicate information from DESCRIPTION section.
 
68
C           (WRB)
 
69
C   920501  Reformatted the REFERENCES section.  (WRB)
 
70
C***END PROLOGUE  DTRSL
 
71
      INTEGER LDT,N,JOB,INFO
 
72
      DOUBLE PRECISION T(LDT,*),B(*)
 
73
C
 
74
C
 
75
      DOUBLE PRECISION DDOT,TEMP
 
76
      INTEGER CASE,J,JJ
 
77
C***FIRST EXECUTABLE STATEMENT  DTRSL
 
78
C
 
79
C        CHECK FOR ZERO DIAGONAL ELEMENTS.
 
80
C
 
81
         DO 10 INFO = 1, N
 
82
            IF (T(INFO,INFO) .EQ. 0.0D0) GO TO 150
 
83
   10    CONTINUE
 
84
         INFO = 0
 
85
C
 
86
C        DETERMINE THE TASK AND GO TO IT.
 
87
C
 
88
         CASE = 1
 
89
         IF (MOD(JOB,10) .NE. 0) CASE = 2
 
90
         IF (MOD(JOB,100)/10 .NE. 0) CASE = CASE + 2
 
91
         GO TO (20,50,80,110), CASE
 
92
C
 
93
C        SOLVE T*X=B FOR T LOWER TRIANGULAR
 
94
C
 
95
   20    CONTINUE
 
96
            B(1) = B(1)/T(1,1)
 
97
            IF (N .LT. 2) GO TO 40
 
98
            DO 30 J = 2, N
 
99
               TEMP = -B(J-1)
 
100
               CALL DAXPY(N-J+1,TEMP,T(J,J-1),1,B(J),1)
 
101
               B(J) = B(J)/T(J,J)
 
102
   30       CONTINUE
 
103
   40       CONTINUE
 
104
         GO TO 140
 
105
C
 
106
C        SOLVE T*X=B FOR T UPPER TRIANGULAR.
 
107
C
 
108
   50    CONTINUE
 
109
            B(N) = B(N)/T(N,N)
 
110
            IF (N .LT. 2) GO TO 70
 
111
            DO 60 JJ = 2, N
 
112
               J = N - JJ + 1
 
113
               TEMP = -B(J+1)
 
114
               CALL DAXPY(J,TEMP,T(1,J+1),1,B(1),1)
 
115
               B(J) = B(J)/T(J,J)
 
116
   60       CONTINUE
 
117
   70       CONTINUE
 
118
         GO TO 140
 
119
C
 
120
C        SOLVE TRANS(T)*X=B FOR T LOWER TRIANGULAR.
 
121
C
 
122
   80    CONTINUE
 
123
            B(N) = B(N)/T(N,N)
 
124
            IF (N .LT. 2) GO TO 100
 
125
            DO 90 JJ = 2, N
 
126
               J = N - JJ + 1
 
127
               B(J) = B(J) - DDOT(JJ-1,T(J+1,J),1,B(J+1),1)
 
128
               B(J) = B(J)/T(J,J)
 
129
   90       CONTINUE
 
130
  100       CONTINUE
 
131
         GO TO 140
 
132
C
 
133
C        SOLVE TRANS(T)*X=B FOR T UPPER TRIANGULAR.
 
134
C
 
135
  110    CONTINUE
 
136
            B(1) = B(1)/T(1,1)
 
137
            IF (N .LT. 2) GO TO 130
 
138
            DO 120 J = 2, N
 
139
               B(J) = B(J) - DDOT(J-1,T(1,J),1,B(1),1)
 
140
               B(J) = B(J)/T(J,J)
 
141
  120       CONTINUE
 
142
  130       CONTINUE
 
143
  140    CONTINUE
 
144
  150 CONTINUE
 
145
      RETURN
 
146
      END