← Back to branch summary

~ubuntu-branches/ubuntu/lucid/pdl/lucid

« back to all changes in this revision

Viewing changes to Lib/Slatec/slatec/sgedi.f

  • Committer: Bazaar Package Importer
  • Author(s): Ben Gertzfield
  • Date: 2002-04-08 18:47:16 UTC
  • Revision ID: james.westby@ubuntu.com-20020408184716-0hf64dc96kin3htp
Tags: upstream-2.3.2
Import upstream version 2.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/Slatec/slatec/sgedi.f
 
1
*DECK SGEDI
 
2
      SUBROUTINE SGEDI (A, LDA, N, IPVT, DET, WORK, JOB)
 
3
C***BEGIN PROLOGUE  SGEDI
 
4
C***PURPOSE  Compute the determinant and inverse of a matrix using the
 
5
C            factors computed by SGECO or SGEFA.
 
6
C***LIBRARY   SLATEC (LINPACK)
 
7
C***CATEGORY  D2A1, D3A1
 
8
C***TYPE      SINGLE PRECISION (SGEDI-S, DGEDI-D, CGEDI-C)
 
9
C***KEYWORDS  DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX
 
10
C***AUTHOR  Moler, C. B., (U. of New Mexico)
 
11
C***DESCRIPTION
 
12
C
 
13
C     SGEDI computes the determinant and inverse of a matrix
 
14
C     using the factors computed by SGECO or SGEFA.
 
15
C
 
16
C     On Entry
 
17
C
 
18
C        A       REAL(LDA, N)
 
19
C                the output from SGECO or SGEFA.
 
20
C
 
21
C        LDA     INTEGER
 
22
C                the leading dimension of the array  A .
 
23
C
 
24
C        N       INTEGER
 
25
C                the order of the matrix  A .
 
26
C
 
27
C        IPVT    INTEGER(N)
 
28
C                the pivot vector from SGECO or SGEFA.
 
29
C
 
30
C        WORK    REAL(N)
 
31
C                work vector.  Contents destroyed.
 
32
C
 
33
C        JOB     INTEGER
 
34
C                = 11   both determinant and inverse.
 
35
C                = 01   inverse only.
 
36
C                = 10   determinant only.
 
37
C
 
38
C     On Return
 
39
C
 
40
C        A       inverse of original matrix if requested.
 
41
C                Otherwise unchanged.
 
42
C
 
43
C        DET     REAL(2)
 
44
C                determinant of original matrix if requested.
 
45
C                Otherwise not referenced.
 
46
C                Determinant = DET(1) * 10.0**DET(2)
 
47
C                with  1.0 .LE. ABS(DET(1)) .LT. 10.0
 
48
C                or  DET(1) .EQ. 0.0 .
 
49
C
 
50
C     Error Condition
 
51
C
 
52
C        A division by zero will occur if the input factor contains
 
53
C        a zero on the diagonal and the inverse is requested.
 
54
C        It will not occur if the subroutines are called correctly
 
55
C        and if SGECO has set RCOND .GT. 0.0 or SGEFA has set
 
56
C        INFO .EQ. 0 .
 
57
C
 
58
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
59
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
60
C***ROUTINES CALLED  SAXPY, SSCAL, SSWAP
 
61
C***REVISION HISTORY  (YYMMDD)
 
62
C   780814  DATE WRITTEN
 
63
C   890831  Modified array declarations.  (WRB)
 
64
C   890831  REVISION DATE from Version 3.2
 
65
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
66
C   900326  Removed duplicate information from DESCRIPTION section.
 
67
C           (WRB)
 
68
C   920501  Reformatted the REFERENCES section.  (WRB)
 
69
C***END PROLOGUE  SGEDI
 
70
      INTEGER LDA,N,IPVT(*),JOB
 
71
      REAL A(LDA,*),DET(2),WORK(*)
 
72
C
 
73
      REAL T
 
74
      REAL TEN
 
75
      INTEGER I,J,K,KB,KP1,L,NM1
 
76
C***FIRST EXECUTABLE STATEMENT  SGEDI
 
77
C
 
78
C     COMPUTE DETERMINANT
 
79
C
 
80
      IF (JOB/10 .EQ. 0) GO TO 70
 
81
         DET(1) = 1.0E0
 
82
         DET(2) = 0.0E0
 
83
         TEN = 10.0E0
 
84
         DO 50 I = 1, N
 
85
            IF (IPVT(I) .NE. I) DET(1) = -DET(1)
 
86
            DET(1) = A(I,I)*DET(1)
 
87
            IF (DET(1) .EQ. 0.0E0) GO TO 60
 
88
   10       IF (ABS(DET(1)) .GE. 1.0E0) GO TO 20
 
89
               DET(1) = TEN*DET(1)
 
90
               DET(2) = DET(2) - 1.0E0
 
91
            GO TO 10
 
92
   20       CONTINUE
 
93
   30       IF (ABS(DET(1)) .LT. TEN) GO TO 40
 
94
               DET(1) = DET(1)/TEN
 
95
               DET(2) = DET(2) + 1.0E0
 
96
            GO TO 30
 
97
   40       CONTINUE
 
98
   50    CONTINUE
 
99
   60    CONTINUE
 
100
   70 CONTINUE
 
101
C
 
102
C     COMPUTE INVERSE(U)
 
103
C
 
104
      IF (MOD(JOB,10) .EQ. 0) GO TO 150
 
105
         DO 100 K = 1, N
 
106
            A(K,K) = 1.0E0/A(K,K)
 
107
            T = -A(K,K)
 
108
            CALL SSCAL(K-1,T,A(1,K),1)
 
109
            KP1 = K + 1
 
110
            IF (N .LT. KP1) GO TO 90
 
111
            DO 80 J = KP1, N
 
112
               T = A(K,J)
 
113
               A(K,J) = 0.0E0
 
114
               CALL SAXPY(K,T,A(1,K),1,A(1,J),1)
 
115
   80       CONTINUE
 
116
   90       CONTINUE
 
117
  100    CONTINUE
 
118
C
 
119
C        FORM INVERSE(U)*INVERSE(L)
 
120
C
 
121
         NM1 = N - 1
 
122
         IF (NM1 .LT. 1) GO TO 140
 
123
         DO 130 KB = 1, NM1
 
124
            K = N - KB
 
125
            KP1 = K + 1
 
126
            DO 110 I = KP1, N
 
127
               WORK(I) = A(I,K)
 
128
               A(I,K) = 0.0E0
 
129
  110       CONTINUE
 
130
            DO 120 J = KP1, N
 
131
               T = WORK(J)
 
132
               CALL SAXPY(N,T,A(1,J),1,A(1,K),1)
 
133
  120       CONTINUE
 
134
            L = IPVT(K)
 
135
            IF (L .NE. K) CALL SSWAP(N,A(1,K),1,A(1,L),1)
 
136
  130    CONTINUE
 
137
  140    CONTINUE
 
138
  150 CONTINUE
 
139
      RETURN
 
140
      END