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

Annotation of /trunk/linpack/dgefa.f

Parent Directory Parent Directory | Revision Log Revision Log


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