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	aw0a	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