Parent Directory | Revision Log

Revision **1** -
(**show annotations**)
(**download**)

*Fri Oct 29 20:54:12 2004 UTC*
(15 years, 5 months ago)
by *aw0a*

File size: 3477 byte(s)

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 |