trunk/linpack/dgbfa.f

C  dgbfa.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 dgbfa(abd,lda,n,ml,mu,ipvt,info)
      integer lda,n,ml,mu,ipvt(1),info
      double precision abd(lda,1)
c
c     dgbfa factors a double precision band matrix by elimination.
c
c     dgbfa is usually called by dgbco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c
c     on entry
c
c        abd     double precision(lda, n)
c                contains the matrix in band storage.  the columns
c                of the matrix are stored in the columns of  abd  and
c                the diagonals of the matrix are stored in rows
c                ml+1 through 2*ml+mu+1 of  abd .
c                see the comments below for details.
c
c        lda     integer
c                the leading dimension of the array  abd .
c                lda must be .ge. 2*ml + mu + 1 .
c
c        n       integer
c                the order of the original matrix.
c
c        ml      integer
c                number of diagonals below the main diagonal.
c                0 .le. ml .lt. n .
c
c        mu      integer
c                number of diagonals above the main diagonal.
c                0 .le. mu .lt. n .
c                more efficient if  ml .le. mu .
c     on return
c
c        abd     an upper triangular matrix in band storage and
c                the multipliers which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgbsl will divide by zero if
c                     called.  use  rcond  in dgbco for a reliable
c                     indication of singularity.
c
c     band storage
c
c           if  a  is a band matrix, the following program segment
c           will set up the input.
c
c                   ml = (band width below the diagonal)
c                   mu = (band width above the diagonal)
c                   m = ml + mu + 1
c                   do 20 j = 1, n
c                      i1 = max0(1, j-mu)
c                      i2 = min0(n, j+ml)
c                      do 10 i = i1, i2
c                         k = i - j + m
c                         abd(k,j) = a(i,j)
c                10    continue
c                20 continue
c
c           this uses rows  ml+1  through  2*ml+mu+1  of  abd .
c           in addition, the first  ml  rows in  abd  are used for
c           elements generated during the triangularization.
c           the total number of rows needed in  abd  is  2*ml+mu+1 .
c           the  ml+mu by ml+mu  upper left triangle and the
c           ml by ml  lower right triangle are not referenced.
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,dscal,idamax
c     fortran max0,min0
c
c     internal variables
c
      double precision t
      integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
c
c
      m = ml + mu + 1
      info = 0
c
c     zero initial fill-in columns
c
      j0 = mu + 2
      j1 = min0(n,m) - 1
      if (j1 .lt. j0) go to 30
      do 20 jz = j0, j1
         i0 = m + 1 - jz
         do 10 i = i0, ml
            abd(i,jz) = 0.0d0
   10    continue
   20 continue
   30 continue
      jz = j1
      ju = 0
c
c     gaussian elimination with partial pivoting
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 130
      do 120 k = 1, nm1
         kp1 = k + 1
c
c        zero next fill-in column
c
         jz = jz + 1
         if (jz .gt. n) go to 50
         if (ml .lt. 1) go to 50
            do 40 i = 1, ml
               abd(i,jz) = 0.0d0
   40       continue
   50    continue
c
c        find l = pivot index
c
         lm = min0(ml,n-k)
         l = idamax(lm+1,abd(m,k),1) + m - 1
         ipvt(k) = l + k - m
c
c        zero pivot implies this column already triangularized
c
         if (abd(l,k) .eq. 0.0d0) go to 100
c
c           interchange if necessary
c
            if (l .eq. m) go to 60
               t = abd(l,k)
               abd(l,k) = abd(m,k)
               abd(m,k) = t
   60       continue
c
c           compute multipliers
c
            t = -1.0d0/abd(m,k)
            call dscal(lm,t,abd(m+1,k),1)
c
c           row elimination with column indexing
c
            ju = min0(max0(ju,mu+ipvt(k)),n)
            mm = m
            if (ju .lt. kp1) go to 90
            do 80 j = kp1, ju
               l = l - 1
               mm = mm - 1
               t = abd(l,j)
               if (l .eq. mm) go to 70
                  abd(l,j) = abd(mm,j)
                  abd(mm,j) = t
   70          continue
               call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
   80       continue
   90       continue
         go to 110
  100    continue
            info = k
  110    continue
  120 continue
  130 continue
      ipvt(n) = n
      if (abd(m,n) .eq. 0.0d0) info = n
      return
      end
1	C dgbfa.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 dgbfa(abd,lda,n,ml,mu,ipvt,info)
7	integer lda,n,ml,mu,ipvt(1),info
8	double precision abd(lda,1)
9	c
10	c dgbfa factors a double precision band matrix by elimination.
11	c
12	c dgbfa is usually called by dgbco, but it can be called
13	c directly with a saving in time if rcond is not needed.
14	c
15	c on entry
16	c
17	c abd double precision(lda, n)
18	c contains the matrix in band storage. the columns
19	c of the matrix are stored in the columns of abd and
20	c the diagonals of the matrix are stored in rows
21	c ml+1 through 2*ml+mu+1 of abd .
22	c see the comments below for details.
23	c
24	c lda integer
25	c the leading dimension of the array abd .
26	c lda must be .ge. 2*ml + mu + 1 .
27	c
28	c n integer
29	c the order of the original matrix.
30	c
31	c ml integer
32	c number of diagonals below the main diagonal.
33	c 0 .le. ml .lt. n .
34	c
35	c mu integer
36	c number of diagonals above the main diagonal.
37	c 0 .le. mu .lt. n .
38	c more efficient if ml .le. mu .
39	c on return
40	c
41	c abd an upper triangular matrix in band storage and
42	c the multipliers which were used to obtain it.
43	c the factorization can be written a = l*u where
44	c l is a product of permutation and unit lower
45	c triangular matrices and u is upper triangular.
46	c
47	c ipvt integer(n)
48	c an integer vector of pivot indices.
49	c
50	c info integer
51	c = 0 normal value.
52	c = k if u(k,k) .eq. 0.0 . this is not an error
53	c condition for this subroutine, but it does
54	c indicate that dgbsl will divide by zero if
55	c called. use rcond in dgbco for a reliable
56	c indication of singularity.
57	c
58	c band storage
59	c
60	c if a is a band matrix, the following program segment
61	c will set up the input.
62	c
63	c ml = (band width below the diagonal)
64	c mu = (band width above the diagonal)
65	c m = ml + mu + 1
66	c do 20 j = 1, n
67	c i1 = max0(1, j-mu)
68	c i2 = min0(n, j+ml)
69	c do 10 i = i1, i2
70	c k = i - j + m
71	c abd(k,j) = a(i,j)
72	c 10 continue
73	c 20 continue
74	c
75	c this uses rows ml+1 through 2*ml+mu+1 of abd .
76	c in addition, the first ml rows in abd are used for
77	c elements generated during the triangularization.
78	c the total number of rows needed in abd is 2*ml+mu+1 .
79	c the ml+mu by ml+mu upper left triangle and the
80	c ml by ml lower right triangle are not referenced.
81	c
82	c linpack. this version dated 08/14/78 .
83	c cleve moler, university of new mexico, argonne national lab.
84	c
85	c subroutines and functions
86	c
87	c blas daxpy,dscal,idamax
88	c fortran max0,min0
89	c
90	c internal variables
91	c
92	double precision t
93	integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
94	c
95	c
96	m = ml + mu + 1
97	info = 0
98	c
99	c zero initial fill-in columns
100	c
101	j0 = mu + 2
102	j1 = min0(n,m) - 1
103	if (j1 .lt. j0) go to 30
104	do 20 jz = j0, j1
105	i0 = m + 1 - jz
106	do 10 i = i0, ml
107	abd(i,jz) = 0.0d0
108	10 continue
109	20 continue
110	30 continue
111	jz = j1
112	ju = 0
113	c
114	c gaussian elimination with partial pivoting
115	c
116	nm1 = n - 1
117	if (nm1 .lt. 1) go to 130
118	do 120 k = 1, nm1
119	kp1 = k + 1
120	c
121	c zero next fill-in column
122	c
123	jz = jz + 1
124	if (jz .gt. n) go to 50
125	if (ml .lt. 1) go to 50
126	do 40 i = 1, ml
127	abd(i,jz) = 0.0d0
128	40 continue
129	50 continue
130	c
131	c find l = pivot index
132	c
133	lm = min0(ml,n-k)
134	l = idamax(lm+1,abd(m,k),1) + m - 1
135	ipvt(k) = l + k - m
136	c
137	c zero pivot implies this column already triangularized
138	c
139	if (abd(l,k) .eq. 0.0d0) go to 100
140	c
141	c interchange if necessary
142	c
143	if (l .eq. m) go to 60
144	t = abd(l,k)
145	abd(l,k) = abd(m,k)
146	abd(m,k) = t
147	60 continue
148	c
149	c compute multipliers
150	c
151	t = -1.0d0/abd(m,k)
152	call dscal(lm,t,abd(m+1,k),1)
153	c
154	c row elimination with column indexing
155	c
156	ju = min0(max0(ju,mu+ipvt(k)),n)
157	mm = m
158	if (ju .lt. kp1) go to 90
159	do 80 j = kp1, ju
160	l = l - 1
161	mm = mm - 1
162	t = abd(l,j)
163	if (l .eq. mm) go to 70
164	abd(l,j) = abd(mm,j)
165	abd(mm,j) = t
166	70 continue
167	call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
168	80 continue
169	90 continue
170	go to 110
171	100 continue
172	info = k
173	110 continue
174	120 continue
175	130 continue
176	ipvt(n) = n
177	if (abd(m,n) .eq. 0.0d0) info = n
178	return
179	end