V3FIT
li_gefa.f
1  SUBROUTINE sgefa1 (A, LDA, N, IPVT, INFO)
2  USE liprec, ONLY: wp => sp
3  IMPLICIT NONE
4 C-----------------------------------------------
5 C D U M M Y A R G U M E N T S
6 C-----------------------------------------------
7  INTEGER LDA, N, INFO
8  INTEGER, DIMENSION(N) :: IPVT
9  REAL(WP), DIMENSION(LDA,N) :: A
10 
11  CALL sgefa (a, lda, n, ipvt, info)
12 
13  END SUBROUTINE sgefa1
14 
15  SUBROUTINE dgefa1 (A, LDA, N, IPVT, INFO)
16  USE liprec, ONLY: wp => dp
17  IMPLICIT NONE
18 C-----------------------------------------------
19 C D U M M Y A R G U M E N T S
20 C-----------------------------------------------
21  INTEGER LDA, N, INFO
22  INTEGER, DIMENSION(N) :: IPVT
23  REAL(WP), DIMENSION(LDA,N) :: A
24 
25  CALL dgefa (a, lda, n, ipvt, info)
26 
27  END SUBROUTINE dgefa1
28 
29 #ifndef CRAY
30  SUBROUTINE sgefa (A, LDA, N, IPVT, INFO)
31  USE liprec, ONLY: wp => sp
32  IMPLICIT NONE
33 C-----------------------------------------------
34 C D U M M Y A R G U M E N T S
35 C-----------------------------------------------
36  INTEGER LDA, N, INFO
37  INTEGER, DIMENSION(N) :: IPVT
38  REAL(WP), DIMENSION(LDA,N) :: A
39 C-----------------------------------------------
40 C L O C A L V A R I A B L E S
41 C-----------------------------------------------
42  REAL(WP), PARAMETER :: ZERO = 0, one = 1
43  INTEGER :: J, K, KP1, L, NM1
44  INTEGER, DIMENSION(1) :: ISAMAX
45  REAL(WP) :: ELEMENT
46 C-----------------------------------------------
47 c
48 c sgefa factors a real matrix by gaussian elimination.
49 c
50 c sgefa is usually called by sgeco, but it can be called
51 c directly with a saving in time if rcond is not needed.
52 c (time for sgeco) = (1 + 9/n)*(time for sgefa) .
53 c
54 c on ENTRY
55 c
56 c a REAL(lda, n)
57 c the matrix to be factored.
58 c
59 c lda INTEGER
60 c the leading dimension of the array a .
61 c
62 c n INTEGER
63 c the order of the matrix a .
64 c
65 c on RETURN
66 c
67 c a an upper triangular matrix and the multipliers
68 c which were used to obtain it.
69 c the factorization can be written a = l*u where
70 c l is a product of permutation and unit lower
71 c triangular matrices and u is upper triangular.
72 c
73 c ipvt INTEGER(n)
74 c an integer vector of pivot indices.
75 c
76 c info integer
77 c = 0 normal value.
78 c = k if u(k,k).eq.0.0 . this is not an error
79 c condition for this subroutine, but it does
80 c indicate that sgesl or sgedi will divide by zero
81 c if called. use rcond in sgeco for a reliable
82 c indication of singularity.
83 c
84 c linpack. this version dated 08/14/78 .
85 c cleve moler, university of new mexico, argonne national lab.
86 c
87 c
88 c
89 c gaussian elimination with partial pivoting
90 c
91  info = 0
92  nm1 = n - 1
93  IF (nm1 .ge. 1) THEN
94  DO k = 1, nm1
95  kp1 = k + 1
96 c
97 c find l = pivot index
98 c
99  isamax = maxloc(abs(a(k:n,k)))
100  l = isamax(1) + k - 1
101 C l = IDAMAX(n-k+1,a(k,k),1) + k - 1
102  ipvt(k) = l
103 c
104 
105 c zero pivot implies this column already triangularized
106 c
107  IF (a(l,k) .ne. zero) THEN
108 c
109 c interchange if necessary
110 c
111  IF (l .ne. k) THEN
112  element = a(l,k)
113  a(l,k) = a(k,k)
114  a(k,k) = element
115  END IF
116 c
117 c compute multipliers
118 c
119  element = -one/a(k,k)
120 ! CALL bla_scal(n-k,element,a(k+1,k),1)
121  a(k+1:n,k) = element*a(k+1:n,k)
122 c
123 c row elimination with column indexing
124 c
125  DO j = kp1, n
126  element = a(l,j)
127  IF (l .ne. k) THEN
128  a(l,j) = a(k,j)
129  a(k,j) = element
130  END IF
131 ! CALL bla_axpy(n-k,element,a(k+1,k),1,a(k+1,j),1)
132  a(k+1:n,j) = a(k+1:n,j) + element*a(k+1:n,k)
133  END DO
134 
135  ELSE
136  info = k
137  END IF
138  END DO
139  END IF
140  ipvt(n) = n
141  IF (a(n,n) .eq. zero) info = n
142 
143  END SUBROUTINE sgefa
144 
145  SUBROUTINE dgefa (A, LDA, N, IPVT, INFO)
146  USE liprec, ONLY: wp => dp
147  IMPLICIT NONE
148 C-----------------------------------------------
149 C D U M M Y A R G U M E N T S
150 C-----------------------------------------------
151  INTEGER LDA, N, INFO
152  INTEGER, DIMENSION(N) :: IPVT
153  REAL(WP), DIMENSION(LDA,N) :: A
154 C-----------------------------------------------
155 C L O C A L V A R I A B L E S
156 C-----------------------------------------------
157  REAL(WP), PARAMETER :: ZERO = 0, one = 1
158  INTEGER :: J, K, KP1, L, NM1
159  INTEGER, DIMENSION(1) :: ISAMAX
160  REAL(WP) :: ELEMENT
161 C-----------------------------------------------
162 c
163 c sgefa factors a REAL matrix by gaussian elimination.
164 c
165 c sgefa is usually called by sgeco, but it can be called
166 c directly with a saving in time IF rcond is not needed.
167 c (time for sgeco) = (1 + 9/n)*(time for sgefa) .
168 c
169 c on ENTRY
170 c
171 c a REAL(lda, n)
172 c the matrix to be factored.
173 c
174 c lda INTEGER
175 c the leading DIMENSION of the array a .
176 c
177 c n INTEGER
178 c the order of the matrix a .
179 c
180 c on RETURN
181 c
182 c a an upper triangular matrix and the multipliers
183 c which were used to obtain it.
184 c the factorization can be written a = l*u WHERE
185 c l is a product of permutation and unit lower
186 c triangular matrices and u is upper triangular.
187 c
188 c ipvt INTEGER(n)
189 c an INTEGER vector of pivot indices.
190 c
191 c info integer
192 c = 0 normal value.
193 c = k if u(k,k).eq.0.0 . this is not an error
194 c condition for this subroutine, but it does
195 c indicate that sgesl or sgedi will divide by zero
196 c if called. use rcond in sgeco for a reliable
197 c indication of singularity.
198 c
199 c linpack. this version dated 08/14/78 .
200 c cleve moler, university of new mexico, argonne national lab.
201 c
202 c
203 c
204 c gaussian elimination with partial pivoting
205 c
206  info = 0
207  nm1 = n - 1
208  IF (nm1 .ge. 1) THEN
209  DO k = 1, nm1
210  kp1 = k + 1
211 c
212 c find l = pivot INDEX
213 c
214  isamax = maxloc(abs(a(k:n,k)))
215  l = isamax(1) + k - 1
216 C l = IDAMAX(n-k+1,a(k,k),1) + k - 1
217  ipvt(k) = l
218 c
219 
220 c zero pivot implies this column already triangularized
221 c
222  IF (a(l,k) .ne. zero) THEN
223 c
224 c interchange if necessary
225 c
226  IF (l .ne. k) THEN
227  element = a(l,k)
228  a(l,k) = a(k,k)
229  a(k,k) = element
230  END IF
231 c
232 c compute multipliers
233 c
234  element = -one/a(k,k)
235 ! CALL bla_scal(n-k,element,a(k+1,k),1)
236  a(k+1:n,k) = element*a(k+1:n,k)
237 c
238 c row elimination with column indexing
239 c
240  DO j = kp1, n
241  element = a(l,j)
242  IF (l .ne. k) THEN
243  a(l,j) = a(k,j)
244  a(k,j) = element
245  END IF
246 ! CALL bla_axpy(n-k,element,a(k+1,k),1,a(k+1,j),1)
247  a(k+1:n,j) = a(k+1:n,j) + element*a(k+1:n,k)
248  END DO
249 
250  ELSE
251  info = k
252  END IF
253  END DO
254  END IF
255  ipvt(n) = n
256  IF (a(n,n) .eq. zero) info = n
257 
258  END SUBROUTINE dgefa
259 #endif