trunk/linpack/dgbsl.f

C  dgbsl.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 dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
      integer lda,n,ml,mu,ipvt(1),job
      double precision abd(lda,1),b(1)
c
c     dgbsl solves the double precision band system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgbco or dgbfa.
c
c     on entry
c
c        abd     double precision(lda, n)
c                the output from dgbco or dgbfa.
c
c        lda     integer
c                the leading dimension of the array  abd .
c
c        n       integer
c                the order of the original matrix.
c
c        ml      integer
c                number of diagonals below the main diagonal.
c
c        mu      integer
c                number of diagonals above the main diagonal.
c
c        ipvt    integer(n)
c                the pivot vector from dgbco or dgbfa.
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 dgbco has set rcond .gt. 0.0
c        or dgbfa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgbsl(abd,lda,n,ml,mu,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     fortran min0
c
c     internal variables
c
      double precision ddot,t
      integer k,kb,l,la,lb,lm,m,nm1
c
      m = mu + ml + 1
      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 (ml .eq. 0) go to 30
         if (nm1 .lt. 1) go to 30
            do 20 k = 1, nm1
               lm = min0(ml,n-k)
               l = ipvt(k)
               t = b(l)
               if (l .eq. k) go to 10
                  b(l) = b(k)
                  b(k) = t
   10          continue
               call daxpy(lm,t,abd(m+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)/abd(m,k)
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = -b(k)
            call daxpy(lm,t,abd(la,k),1,b(lb),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
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = ddot(lm,abd(la,k),1,b(lb),1)
            b(k) = (b(k) - t)/abd(m,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (ml .eq. 0) go to 90
         if (nm1 .lt. 1) go to 90
            do 80 kb = 1, nm1
               k = n - kb
               lm = min0(ml,n-k)
               b(k) = b(k) + ddot(lm,abd(m+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	C dgbsl.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 dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
7	integer lda,n,ml,mu,ipvt(1),job
8	double precision abd(lda,1),b(1)
9	c
10	c dgbsl solves the double precision band system
11	c a * x = b or trans(a) * x = b
12	c using the factors computed by dgbco or dgbfa.
13	c
14	c on entry
15	c
16	c abd double precision(lda, n)
17	c the output from dgbco or dgbfa.
18	c
19	c lda integer
20	c the leading dimension of the array abd .
21	c
22	c n integer
23	c the order of the original matrix.
24	c
25	c ml integer
26	c number of diagonals below the main diagonal.
27	c
28	c mu integer
29	c number of diagonals above the main diagonal.
30	c
31	c ipvt integer(n)
32	c the pivot vector from dgbco or dgbfa.
33	c
34	c b double precision(n)
35	c the right hand side vector.
36	c
37	c job integer
38	c = 0 to solve a*x = b ,
39	c = nonzero to solve trans(a)*x = b , where
40	c trans(a) is the transpose.
41	c
42	c on return
43	c
44	c b the solution vector x .
45	c
46	c error condition
47	c
48	c a division by zero will occur if the input factor contains a
49	c zero on the diagonal. technically this indicates singularity
50	c but it is often caused by improper arguments or improper
51	c setting of lda . it will not occur if the subroutines are
52	c called correctly and if dgbco has set rcond .gt. 0.0
53	c or dgbfa has set info .eq. 0 .
54	c
55	c to compute inverse(a) * c where c is a matrix
56	c with p columns
57	c call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
58	c if (rcond is too small) go to ...
59	c do 10 j = 1, p
60	c call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
61	c 10 continue
62	c
63	c linpack. this version dated 08/14/78 .
64	c cleve moler, university of new mexico, argonne national lab.
65	c
66	c subroutines and functions
67	c
68	c blas daxpy,ddot
69	c fortran min0
70	c
71	c internal variables
72	c
73	double precision ddot,t
74	integer k,kb,l,la,lb,lm,m,nm1
75	c
76	m = mu + ml + 1
77	nm1 = n - 1
78	if (job .ne. 0) go to 50
79	c
80	c job = 0 , solve a * x = b
81	c first solve l*y = b
82	c
83	if (ml .eq. 0) go to 30
84	if (nm1 .lt. 1) go to 30
85	do 20 k = 1, nm1
86	lm = min0(ml,n-k)
87	l = ipvt(k)
88	t = b(l)
89	if (l .eq. k) go to 10
90	b(l) = b(k)
91	b(k) = t
92	10 continue
93	call daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
94	20 continue
95	30 continue
96	c
97	c now solve u*x = y
98	c
99	do 40 kb = 1, n
100	k = n + 1 - kb
101	b(k) = b(k)/abd(m,k)
102	lm = min0(k,m) - 1
103	la = m - lm
104	lb = k - lm
105	t = -b(k)
106	call daxpy(lm,t,abd(la,k),1,b(lb),1)
107	40 continue
108	go to 100
109	50 continue
110	c
111	c job = nonzero, solve trans(a) * x = b
112	c first solve trans(u)*y = b
113	c
114	do 60 k = 1, n
115	lm = min0(k,m) - 1
116	la = m - lm
117	lb = k - lm
118	t = ddot(lm,abd(la,k),1,b(lb),1)
119	b(k) = (b(k) - t)/abd(m,k)
120	60 continue
121	c
122	c now solve trans(l)*x = y
123	c
124	if (ml .eq. 0) go to 90
125	if (nm1 .lt. 1) go to 90
126	do 80 kb = 1, nm1
127	k = n - kb
128	lm = min0(ml,n-k)
129	b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
130	l = ipvt(k)
131	if (l .eq. k) go to 70
132	t = b(l)
133	b(l) = b(k)
134	b(k) = t
135	70 continue
136	80 continue
137	90 continue
138	100 continue
139	return
140	end