trunk/linpack/dgesl.f

C  dgesl.f
C  is freely available from netlib. It is not subject to any GNU License
C  set by the authors of the ASCEND math programming system.
C  $Date: 1996/04/30 18:17:11 $ $Revision: 1.1.1.1 $
C
      subroutine dgesl(a,lda,n,ipvt,b,job)
      integer lda,n,ipvt(1),job
      double precision a(lda,1),b(1)
c
c     dgesl solves the double precision system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgeco or dgefa.
c
c     on entry
c
c        a       double precision(lda, n)
c                the output from dgeco or dgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from dgeco or dgefa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgeco has set rcond .gt. 0.0
c        or dgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c
c     internal variables
c
      double precision ddot,t
      integer k,kb,l,nm1
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call daxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            t = ddot(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end
1	aw0a	1	C dgesl.f
2			C is freely available from netlib. It is not subject to any GNU License
3			C set by the authors of the ASCEND math programming system.
4			C $Date: 1996/04/30 18:17:11 $ $Revision: 1.1.1.1 $
5			C
6			subroutine dgesl(a,lda,n,ipvt,b,job)
7			integer lda,n,ipvt(1),job
8			double precision a(lda,1),b(1)
9			c
10			c dgesl solves the double precision system
11			c a * x = b or trans(a) * x = b
12			c using the factors computed by dgeco or dgefa.
13			c
14			c on entry
15			c
16			c a double precision(lda, n)
17			c the output from dgeco or dgefa.
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 dgeco or dgefa.
27			c
28			c b double precision(n)
29			c the right hand side vector.
30			c
31			c job integer
32			c = 0 to solve a*x = b ,
33			c = nonzero to solve trans(a)*x = b where
34			c trans(a) is the transpose.
35			c
36			c on return
37			c
38			c b the solution vector x .
39			c
40			c error condition
41			c
42			c a division by zero will occur if the input factor contains a
43			c zero on the diagonal. technically this indicates singularity
44			c but it is often caused by improper arguments or improper
45			c setting of lda . it will not occur if the subroutines are
46			c called correctly and if dgeco has set rcond .gt. 0.0
47			c or dgefa has set info .eq. 0 .
48			c
49			c to compute inverse(a) * c where c is a matrix
50			c with p columns
51			c call dgeco(a,lda,n,ipvt,rcond,z)
52			c if (rcond is too small) go to ...
53			c do 10 j = 1, p
54			c call dgesl(a,lda,n,ipvt,c(1,j),0)
55			c 10 continue
56			c
57			c linpack. this version dated 08/14/78 .
58			c cleve moler, university of new mexico, argonne national lab.
59			c
60			c subroutines and functions
61			c
62			c blas daxpy,ddot
63			c
64			c internal variables
65			c
66			double precision ddot,t
67			integer k,kb,l,nm1
68			c
69			nm1 = n - 1
70			if (job .ne. 0) go to 50
71			c
72			c job = 0 , solve a * x = b
73			c first solve l*y = b
74			c
75			if (nm1 .lt. 1) go to 30
76			do 20 k = 1, nm1
77			l = ipvt(k)
78			t = b(l)
79			if (l .eq. k) go to 10
80			b(l) = b(k)
81			b(k) = t
82			10 continue
83			call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
84			20 continue
85			30 continue
86			c
87			c now solve u*x = y
88			c
89			do 40 kb = 1, n
90			k = n + 1 - kb
91			b(k) = b(k)/a(k,k)
92			t = -b(k)
93			call daxpy(k-1,t,a(1,k),1,b(1),1)
94			40 continue
95			go to 100
96			50 continue
97			c
98			c job = nonzero, solve trans(a) * x = b
99			c first solve trans(u)*y = b
100			c
101			do 60 k = 1, n
102			t = ddot(k-1,a(1,k),1,b(1),1)
103			b(k) = (b(k) - t)/a(k,k)
104			60 continue
105			c
106			c now solve trans(l)*x = y
107			c
108			if (nm1 .lt. 1) go to 90
109			do 80 kb = 1, nm1
110			k = n - kb
111			b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
112			l = ipvt(k)
113			if (l .eq. k) go to 70
114			t = b(l)
115			b(l) = b(k)
116			b(k) = t
117			70 continue
118			80 continue
119			90 continue
120			100 continue
121			return
122			end