/[ascend]/trunk/linpack/dgefa.f
ViewVC logotype

Contents of /trunk/linpack/dgefa.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (16 years, 1 month ago) by aw0a
File size: 3477 byte(s)
Setting up web subdirectory in repository
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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22