← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/control/wgedi.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/control/wgedi.f
 
1
      subroutine wgedi(ar,ai,lda,n,ipvt,detr,deti,workr,worki,job)
 
2
c     Copyright INRIA
 
3
      integer lda,n,ipvt(*),job
 
4
      double precision ar(lda,*),ai(lda,*),detr(2),deti(2),workr(*),
 
5
     *                 worki(*)
 
6
c!purpose
 
7
c
 
8
c     wgedi computes the determinant and inverse of a matrix
 
9
c     using the factors computed by wgeco or wgefa.
 
10
c
 
11
c!calling sequence
 
12
c
 
13
c      subroutine wgedi(ar,ai,lda,n,ipvt,detr,deti,workr,worki,job)
 
14
c     on entry
 
15
c
 
16
c        a       double-complex(lda, n)
 
17
c                the output from wgeco or wgefa.
 
18
c
 
19
c        lda     integer
 
20
c                the leading dimension of the array  a .
 
21
c
 
22
c        n       integer
 
23
c                the order of the matrix  a .
 
24
c
 
25
c        ipvt    integer(n)
 
26
c                the pivot vector from wgeco or wgefa.
 
27
c
 
28
c        work    double-complex(n)
 
29
c                work vector.  contents destroyed.
 
30
c
 
31
c        job     integer
 
32
c                = 11   both determinant and inverse.
 
33
c                = 01   inverse only.
 
34
c                = 10   determinant only.
 
35
c
 
36
c     on return
 
37
c
 
38
c        a       inverse of original matrix if requested.
 
39
c                otherwise unchanged.
 
40
c
 
41
c        det     double-complex(2)
 
42
c                determinant of original matrix if requested.
 
43
c                otherwise not referenced.
 
44
c                determinant = det(1) * 10.0**det(2)
 
45
c                with  1.0 .le. cabs1(det(1) .lt. 10.0
 
46
c                or  det(1) .eq. 0.0 .
 
47
c
 
48
c     error condition
 
49
c
 
50
c        a division by zero will occur if the input factor contains
 
51
c        a zero on the diagonal and the inverse is requested.
 
52
c        it will not occur if the subroutines are called correctly
 
53
c        and if wgeco has set rcond .gt. 0.0 or wgefa has set
 
54
c        info .eq. 0 .
 
55
c
 
56
c!originator
 
57
c     linpack. this version dated 07/01/79 .
 
58
c     cleve moler, university of new mexico, argonne national lab.
 
59
c
 
60
c!auxiliary routines
 
61
c
 
62
c     blas waxpy,wscal,wswap
 
63
c     fortran abs,mod
 
64
c
 
65
c!
 
66
c     internal variables
 
67
c
 
68
      double precision tr,ti
 
69
      double precision ten
 
70
      integer i,j,k,kb,kp1,l,nm1
 
71
c
 
72
      double precision zdumr,zdumi
 
73
      double precision cabs1
 
74
      cabs1(zdumr,zdumi) = abs(zdumr) + abs(zdumi)
 
75
c
 
76
c     compute determinant
 
77
c
 
78
      if (job/10 .eq. 0) go to 80
 
79
         detr(1) = 1.0d+0
 
80
         deti(1) = 0.0d+0
 
81
         detr(2) = 0.0d+0
 
82
         deti(2) = 0.0d+0
 
83
         ten = 10.0d+0
 
84
         do 60 i = 1, n
 
85
            if (ipvt(i) .eq. i) go to 10
 
86
               detr(1) = -detr(1)
 
87
               deti(1) = -deti(1)
 
88
   10       continue
 
89
            call wmul(ar(i,i),ai(i,i),detr(1),deti(1),detr(1),deti(1))
 
90
c           ...exit
 
91
c        ...exit
 
92
            if (cabs1(detr(1),deti(1)) .eq. 0.0d+0) go to 70
 
93
   20       if (cabs1(detr(1),deti(1)) .ge. 1.0d+0) go to 30
 
94
               detr(1) = ten*detr(1)
 
95
               deti(1) = ten*deti(1)
 
96
               detr(2) = detr(2) - 1.0d+0
 
97
               deti(2) = deti(2) - 0.0d+0
 
98
            go to 20
 
99
   30       continue
 
100
   40       if (cabs1(detr(1),deti(1)) .lt. ten) go to 50
 
101
               detr(1) = detr(1)/ten
 
102
               deti(1) = deti(1)/ten
 
103
               detr(2) = detr(2) + 1.0d+0
 
104
               deti(2) = deti(2) + 0.0d+0
 
105
            go to 40
 
106
   50       continue
 
107
   60    continue
 
108
   70    continue
 
109
   80 continue
 
110
c
 
111
c     compute inverse(u)
 
112
c
 
113
      if (mod(job,10) .eq. 0) go to 160
 
114
         do 110 k = 1, n
 
115
            call wdiv(1.0d+0,0.0d+0,ar(k,k),ai(k,k),ar(k,k),ai(k,k))
 
116
            tr = -ar(k,k)
 
117
            ti = -ai(k,k)
 
118
            call wscal(k-1,tr,ti,ar(1,k),ai(1,k),1)
 
119
            kp1 = k + 1
 
120
            if (n .lt. kp1) go to 100
 
121
            do 90 j = kp1, n
 
122
               tr = ar(k,j)
 
123
               ti = ai(k,j)
 
124
               ar(k,j) = 0.0d+0
 
125
               ai(k,j) = 0.0d+0
 
126
               call waxpy(k,tr,ti,ar(1,k),ai(1,k),1,ar(1,j),ai(1,j),1)
 
127
   90       continue
 
128
  100       continue
 
129
  110    continue
 
130
c
 
131
c        form inverse(u)*inverse(l)
 
132
c
 
133
         nm1 = n - 1
 
134
         if (nm1 .lt. 1) go to 150
 
135
         do 140 kb = 1, nm1
 
136
            k = n - kb
 
137
            kp1 = k + 1
 
138
            do 120 i = kp1, n
 
139
               workr(i) = ar(i,k)
 
140
               worki(i) = ai(i,k)
 
141
               ar(i,k) = 0.0d+0
 
142
               ai(i,k) = 0.0d+0
 
143
  120       continue
 
144
            do 130 j = kp1, n
 
145
               tr = workr(j)
 
146
               ti = worki(j)
 
147
               call waxpy(n,tr,ti,ar(1,j),ai(1,j),1,ar(1,k),ai(1,k),1)
 
148
  130       continue
 
149
            l = ipvt(k)
 
150
            if (l .ne. k)
 
151
     *         call wswap(n,ar(1,k),ai(1,k),1,ar(1,l),ai(1,l),1)
 
152
  140    continue
 
153
  150    continue
 
154
  160 continue
 
155
      return
 
156
      end