Parent Directory | 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)

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 |