trunk/linpack/dgefa.f

C  dgefa.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: 1998/07/06 21:46:42 $ $Revision: 1.3 $
C
      subroutine dgefa(a,lda,n,ipvt,info)
      integer lda,n,ipvt(1),info
      double precision a(lda,1)
c
c     dgefa factors a double precision matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the matrix to be factored.
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     on return
c
c        a       an upper triangular matrix and the multipliers
c                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 dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
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
c     internal variables
c
      double precision t
      integer idamax,j,k,kp1,l,nm1
c
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0d0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -1.0d0/a(k,k)
            call dscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0d0) info = n
      return
      end

C
C The following is a routine to debug the input to
C dgefa from lsode. It should not be called in production code.
C
      subroutine ascmdump(a,lda,n,ipvt,info)
      integer lda,n,ipvt(1),info
      double precision a(lda,1)
      integer j,k

      open(unit=22,file='/tmp/fjex.m',status='UNKNOWN')
      do 1000 j=1,n
      do 1100 k=1,n
        write(22,*) j, k, a(j,k)
1100  continue
1000  continue
      close(22)
      return
      END
1	C dgefa.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: 1998/07/06 21:46:42 $ $Revision: 1.3 $
5	C
6	subroutine dgefa(a,lda,n,ipvt,info)
7	integer lda,n,ipvt(1),info
8	double precision a(lda,1)
9	c
10	c dgefa factors a double precision matrix by gaussian elimination.
11	c
12	c dgefa is usually called by dgeco, but it can be called
13	c directly with a saving in time if rcond is not needed.
14	c (time for dgeco) = (1 + 9/n)*(time for dgefa) .
15	c
16	c on entry
17	c
18	c a double precision(lda, n)
19	c the matrix to be factored.
20	c
21	c lda integer
22	c the leading dimension of the array a .
23	c
24	c n integer
25	c the order of the matrix a .
26	c
27	c on return
28	c
29	c a an upper triangular matrix and the multipliers
30	c which were used to obtain it.
31	c the factorization can be written a = l*u where
32	c l is a product of permutation and unit lower
33	c triangular matrices and u is upper triangular.
34	c
35	c ipvt integer(n)
36	c an integer vector of pivot indices.
37	c
38	c info integer
39	c = 0 normal value.
40	c = k if u(k,k) .eq. 0.0 . this is not an error
41	c condition for this subroutine, but it does
42	c indicate that dgesl or dgedi will divide by zero
43	c if called. use rcond in dgeco for a reliable
44	c indication of singularity.
45	c
46	c linpack. this version dated 08/14/78 .
47	c cleve moler, university of new mexico, argonne national lab.
48	c
49	c subroutines and functions
50	c
51	c blas daxpy,dscal,idamax
52	c
53	c internal variables
54	c
55	double precision t
56	integer idamax,j,k,kp1,l,nm1
57	c
58	c
59	c gaussian elimination with partial pivoting
60	c
61	info = 0
62	nm1 = n - 1
63	if (nm1 .lt. 1) go to 70
64	do 60 k = 1, nm1
65	kp1 = k + 1
66	c
67	c find l = pivot index
68	c
69	l = idamax(n-k+1,a(k,k),1) + k - 1
70	ipvt(k) = l
71	c
72	c zero pivot implies this column already triangularized
73	c
74	if (a(l,k) .eq. 0.0d0) go to 40
75	c
76	c interchange if necessary
77	c
78	if (l .eq. k) go to 10
79	t = a(l,k)
80	a(l,k) = a(k,k)
81	a(k,k) = t
82	10 continue
83	c
84	c compute multipliers
85	c
86	t = -1.0d0/a(k,k)
87	call dscal(n-k,t,a(k+1,k),1)
88	c
89	c row elimination with column indexing
90	c
91	do 30 j = kp1, n
92	t = a(l,j)
93	if (l .eq. k) go to 20
94	a(l,j) = a(k,j)
95	a(k,j) = t
96	20 continue
97	call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
98	30 continue
99	go to 50
100	40 continue
101	info = k
102	50 continue
103	60 continue
104	70 continue
105	ipvt(n) = n
106	if (a(n,n) .eq. 0.0d0) info = n
107	return
108	end
109
110	C
111	C The following is a routine to debug the input to
112	C dgefa from lsode. It should not be called in production code.
113	C
114	subroutine ascmdump(a,lda,n,ipvt,info)
115	integer lda,n,ipvt(1),info
116	double precision a(lda,1)
117	integer j,k
118
119	open(unit=22,file='/tmp/fjex.m',status='UNKNOWN')
120	do 1000 j=1,n
121	do 1100 k=1,n
122	write(22,*) j, k, a(j,k)
123	1100 continue
124	1000 continue
125	close(22)
126	return
127	END