source: CIVL/examples/fortran/nek5000/verification/mxm_pencil.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 11.0 KB
Line 
1 subroutine mxm(a, m, b, k, c, n)
2 real a(m,k), b(k,n), c(m,n)
3 real s11, s12, s13, s14, s21, s22, s23, s24
4 real s31, s32, s33, s34, s41, s42, s43, s44
5
6 mresid = mod(m,4)
7 nresid = mod(n,4)
8 m1 = m - mresid + 1
9 n1 = n - nresid + 1
10
11 do i=1,m-mresid,4
12 do j=1,n-nresid,4
13 s11 = 0.0d0
14 s21 = 0.0d0
15 s31 = 0.0d0
16 s41 = 0.0d0
17 s12 = 0.0d0
18 s22 = 0.0d0
19 s32 = 0.0d0
20 s42 = 0.0d0
21 s13 = 0.0d0
22 s23 = 0.0d0
23 s33 = 0.0d0
24 s43 = 0.0d0
25 s14 = 0.0d0
26 s24 = 0.0d0
27 s34 = 0.0d0
28 s44 = 0.0d0
29 do l=1,k
30 s11 = s11 + a(i,l)*b(l,j)
31 s12 = s12 + a(i,l)*b(l,j+1)
32 s13 = s13 + a(i,l)*b(l,j+2)
33 s14 = s14 + a(i,l)*b(l,j+3)
34
35 s21 = s21 + a(i+1,l)*b(l,j)
36 s22 = s22 + a(i+1,l)*b(l,j+1)
37 s23 = s23 + a(i+1,l)*b(l,j+2)
38 s24 = s24 + a(i+1,l)*b(l,j+3)
39
40 s31 = s31 + a(i+2,l)*b(l,j)
41 s32 = s32 + a(i+2,l)*b(l,j+1)
42 s33 = s33 + a(i+2,l)*b(l,j+2)
43 s34 = s34 + a(i+2,l)*b(l,j+3)
44
45 s41 = s41 + a(i+3,l)*b(l,j)
46 s42 = s42 + a(i+3,l)*b(l,j+1)
47 s43 = s43 + a(i+3,l)*b(l,j+2)
48 s44 = s44 + a(i+3,l)*b(l,j+3)
49 enddo
50 c(i,j) = s11
51 c(i,j+1) = s12
52 c(i,j+2) = s13
53 c(i,j+3) = s14
54
55 c(i+1,j) = s21
56 c(i+2,j) = s31
57 c(i+3,j) = s41
58
59 c(i+1,j+1) = s22
60 c(i+2,j+1) = s32
61 c(i+3,j+1) = s42
62
63 c(i+1,j+2) = s23
64 c(i+2,j+2) = s33
65 c(i+3,j+2) = s43
66
67 c(i+1,j+3) = s24
68 c(i+2,j+3) = s34
69 c(i+3,j+3) = s44
70 enddo
71* Residual when n is not multiple of 4
72 if (nresid .ne. 0) then
73 if (nresid .eq. 1) then
74 s11 = 0.0d0
75 s21 = 0.0d0
76 s31 = 0.0d0
77 s41 = 0.0d0
78 do l=1,k
79 s11 = s11 + a(i,l)*b(l,n)
80 s21 = s21 + a(i+1,l)*b(l,n)
81 s31 = s31 + a(i+2,l)*b(l,n)
82 s41 = s41 + a(i+3,l)*b(l,n)
83 enddo
84 c(i,n) = s11
85 c(i+1,n) = s21
86 c(i+2,n) = s31
87 c(i+3,n) = s41
88 elseif (nresid .eq. 2) then
89 s11 = 0.0d0
90 s21 = 0.0d0
91 s31 = 0.0d0
92 s41 = 0.0d0
93 s12 = 0.0d0
94 s22 = 0.0d0
95 s32 = 0.0d0
96 s42 = 0.0d0
97 do l=1,k
98 s11 = s11 + a(i,l)*b(l,j)
99 s12 = s12 + a(i,l)*b(l,j+1)
100
101 s21 = s21 + a(i+1,l)*b(l,j)
102 s22 = s22 + a(i+1,l)*b(l,j+1)
103
104 s31 = s31 + a(i+2,l)*b(l,j)
105 s32 = s32 + a(i+2,l)*b(l,j+1)
106
107 s41 = s41 + a(i+3,l)*b(l,j)
108 s42 = s42 + a(i+3,l)*b(l,j+1)
109 enddo
110 c(i,j) = s11
111 c(i,j+1) = s12
112
113 c(i+1,j) = s21
114 c(i+2,j) = s31
115 c(i+3,j) = s41
116
117 c(i+1,j+1) = s22
118 c(i+2,j+1) = s32
119 c(i+3,j+1) = s42
120 else
121 s11 = 0.0d0
122 s21 = 0.0d0
123 s31 = 0.0d0
124 s41 = 0.0d0
125 s12 = 0.0d0
126 s22 = 0.0d0
127 s32 = 0.0d0
128 s42 = 0.0d0
129 s13 = 0.0d0
130 s23 = 0.0d0
131 s33 = 0.0d0
132 s43 = 0.0d0
133 do l=1,k
134 s11 = s11 + a(i,l)*b(l,j)
135 s12 = s12 + a(i,l)*b(l,j+1)
136 s13 = s13 + a(i,l)*b(l,j+2)
137
138 s21 = s21 + a(i+1,l)*b(l,j)
139 s22 = s22 + a(i+1,l)*b(l,j+1)
140 s23 = s23 + a(i+1,l)*b(l,j+2)
141
142 s31 = s31 + a(i+2,l)*b(l,j)
143 s32 = s32 + a(i+2,l)*b(l,j+1)
144 s33 = s33 + a(i+2,l)*b(l,j+2)
145
146 s41 = s41 + a(i+3,l)*b(l,j)
147 s42 = s42 + a(i+3,l)*b(l,j+1)
148 s43 = s43 + a(i+3,l)*b(l,j+2)
149 enddo
150 c(i,j) = s11
151 c(i+1,j) = s21
152 c(i+2,j) = s31
153 c(i+3,j) = s41
154 c(i,j+1) = s12
155 c(i+1,j+1) = s22
156 c(i+2,j+1) = s32
157 c(i+3,j+1) = s42
158 c(i,j+2) = s13
159 c(i+1,j+2) = s23
160 c(i+2,j+2) = s33
161 c(i+3,j+2) = s43
162 endif
163 endif
164 enddo
165
166* Residual when m is not multiple of 4
167 if (mresid .eq. 0) then
168 return
169 elseif (mresid .eq. 1) then
170 do j=1,n-nresid,4
171 s11 = 0.0d0
172 s12 = 0.0d0
173 s13 = 0.0d0
174 s14 = 0.0d0
175 do l=1,k
176 s11 = s11 + a(m,l)*b(l,j)
177 s12 = s12 + a(m,l)*b(l,j+1)
178 s13 = s13 + a(m,l)*b(l,j+2)
179 s14 = s14 + a(m,l)*b(l,j+3)
180 enddo
181 c(m,j) = s11
182 c(m,j+1) = s12
183 c(m,j+2) = s13
184 c(m,j+3) = s14
185 enddo
186* mresid is 1, check nresid
187 if (nresid .eq. 0) then
188 return
189 elseif (nresid .eq. 1) then
190 s11 = 0.0d0
191 do l=1,k
192 s11 = s11 + a(m,l)*b(l,n)
193 enddo
194 c(m,n) = s11
195 return
196 elseif (nresid .eq. 2) then
197 s11 = 0.0d0
198 s12 = 0.0d0
199 do l=1,k
200 s11 = s11 + a(m,l)*b(l,n-1)
201 s12 = s12 + a(m,l)*b(l,n)
202 enddo
203 c(m,n-1) = s11
204 c(m,n) = s12
205 return
206 else
207 s11 = 0.0d0
208 s12 = 0.0d0
209 s13 = 0.0d0
210 do l=1,k
211 s11 = s11 + a(m,l)*b(l,n-2)
212 s12 = s12 + a(m,l)*b(l,n-1)
213 s13 = s13 + a(m,l)*b(l,n)
214 enddo
215 c(m,n-2) = s11
216 c(m,n-1) = s12
217 c(m,n) = s13
218 return
219 endif
220 elseif (mresid .eq. 2) then
221 do j=1,n-nresid,4
222 s11 = 0.0d0
223 s12 = 0.0d0
224 s13 = 0.0d0
225 s14 = 0.0d0
226 s21 = 0.0d0
227 s22 = 0.0d0
228 s23 = 0.0d0
229 s24 = 0.0d0
230 do l=1,k
231 s11 = s11 + a(m-1,l)*b(l,j)
232 s12 = s12 + a(m-1,l)*b(l,j+1)
233 s13 = s13 + a(m-1,l)*b(l,j+2)
234 s14 = s14 + a(m-1,l)*b(l,j+3)
235
236 s21 = s21 + a(m,l)*b(l,j)
237 s22 = s22 + a(m,l)*b(l,j+1)
238 s23 = s23 + a(m,l)*b(l,j+2)
239 s24 = s24 + a(m,l)*b(l,j+3)
240 enddo
241 c(m-1,j) = s11
242 c(m-1,j+1) = s12
243 c(m-1,j+2) = s13
244 c(m-1,j+3) = s14
245 c(m,j) = s21
246 c(m,j+1) = s22
247 c(m,j+2) = s23
248 c(m,j+3) = s24
249 enddo
250* mresid is 2, check nresid
251 if (nresid .eq. 0) then
252 return
253 elseif (nresid .eq. 1) then
254 s11 = 0.0d0
255 s21 = 0.0d0
256 do l=1,k
257 s11 = s11 + a(m-1,l)*b(l,n)
258 s21 = s21 + a(m,l)*b(l,n)
259 enddo
260 c(m-1,n) = s11
261 c(m,n) = s21
262 return
263 elseif (nresid .eq. 2) then
264 s11 = 0.0d0
265 s21 = 0.0d0
266 s12 = 0.0d0
267 s22 = 0.0d0
268 do l=1,k
269 s11 = s11 + a(m-1,l)*b(l,n-1)
270 s12 = s12 + a(m-1,l)*b(l,n)
271 s21 = s21 + a(m,l)*b(l,n-1)
272 s22 = s22 + a(m,l)*b(l,n)
273 enddo
274 c(m-1,n-1) = s11
275 c(m-1,n) = s12
276 c(m,n-1) = s21
277 c(m,n) = s22
278 return
279 else
280 s11 = 0.0d0
281 s21 = 0.0d0
282 s12 = 0.0d0
283 s22 = 0.0d0
284 s13 = 0.0d0
285 s23 = 0.0d0
286 do l=1,k
287 s11 = s11 + a(m-1,l)*b(l,n-2)
288 s12 = s12 + a(m-1,l)*b(l,n-1)
289 s13 = s13 + a(m-1,l)*b(l,n)
290 s21 = s21 + a(m,l)*b(l,n-2)
291 s22 = s22 + a(m,l)*b(l,n-1)
292 s23 = s23 + a(m,l)*b(l,n)
293 enddo
294 c(m-1,n-2) = s11
295 c(m-1,n-1) = s12
296 c(m-1,n) = s13
297 c(m,n-2) = s21
298 c(m,n-1) = s22
299 c(m,n) = s23
300 return
301 endif
302 else
303* mresid is 3
304 do j=1,n-nresid,4
305 s11 = 0.0d0
306 s21 = 0.0d0
307 s31 = 0.0d0
308
309 s12 = 0.0d0
310 s22 = 0.0d0
311 s32 = 0.0d0
312
313 s13 = 0.0d0
314 s23 = 0.0d0
315 s33 = 0.0d0
316
317 s14 = 0.0d0
318 s24 = 0.0d0
319 s34 = 0.0d0
320
321 do l=1,k
322 s11 = s11 + a(m-2,l)*b(l,j)
323 s12 = s12 + a(m-2,l)*b(l,j+1)
324 s13 = s13 + a(m-2,l)*b(l,j+2)
325 s14 = s14 + a(m-2,l)*b(l,j+3)
326
327 s21 = s21 + a(m-1,l)*b(l,j)
328 s22 = s22 + a(m-1,l)*b(l,j+1)
329 s23 = s23 + a(m-1,l)*b(l,j+2)
330 s24 = s24 + a(m-1,l)*b(l,j+3)
331
332 s31 = s31 + a(m,l)*b(l,j)
333 s32 = s32 + a(m,l)*b(l,j+1)
334 s33 = s33 + a(m,l)*b(l,j+2)
335 s34 = s34 + a(m,l)*b(l,j+3)
336 enddo
337 c(m-2,j) = s11
338 c(m-2,j+1) = s12
339 c(m-2,j+2) = s13
340 c(m-2,j+3) = s14
341
342 c(m-1,j) = s21
343 c(m-1,j+1) = s22
344 c(m-1,j+2) = s23
345 c(m-1,j+3) = s24
346
347 c(m,j) = s31
348 c(m,j+1) = s32
349 c(m,j+2) = s33
350 c(m,j+3) = s34
351 enddo
352* mresid is 3, check nresid
353 if (nresid .eq. 0) then
354 return
355 elseif (nresid .eq. 1) then
356 s11 = 0.0d0
357 s21 = 0.0d0
358 s31 = 0.0d0
359 do l=1,k
360 s11 = s11 + a(m-2,l)*b(l,n)
361 s21 = s21 + a(m-1,l)*b(l,n)
362 s31 = s31 + a(m,l)*b(l,n)
363 enddo
364 c(m-2,n) = s11
365 c(m-1,n) = s21
366 c(m,n) = s31
367 return
368 elseif (nresid .eq. 2) then
369 s11 = 0.0d0
370 s21 = 0.0d0
371 s31 = 0.0d0
372 s12 = 0.0d0
373 s22 = 0.0d0
374 s32 = 0.0d0
375 do l=1,k
376 s11 = s11 + a(m-2,l)*b(l,n-1)
377 s12 = s12 + a(m-2,l)*b(l,n)
378 s21 = s21 + a(m-1,l)*b(l,n-1)
379 s22 = s22 + a(m-1,l)*b(l,n)
380 s31 = s31 + a(m,l)*b(l,n-1)
381 s32 = s32 + a(m,l)*b(l,n)
382 enddo
383 c(m-2,n-1) = s11
384 c(m-2,n) = s12
385 c(m-1,n-1) = s21
386 c(m-1,n) = s22
387 c(m,n-1) = s31
388 c(m,n) = s32
389 return
390 else
391 s11 = 0.0d0
392 s21 = 0.0d0
393 s31 = 0.0d0
394 s12 = 0.0d0
395 s22 = 0.0d0
396 s32 = 0.0d0
397 s13 = 0.0d0
398 s23 = 0.0d0
399 s33 = 0.0d0
400 do l=1,k
401 s11 = s11 + a(m-2,l)*b(l,n-2)
402 s12 = s12 + a(m-2,l)*b(l,n-1)
403 s13 = s13 + a(m-2,l)*b(l,n)
404 s21 = s21 + a(m-1,l)*b(l,n-2)
405 s22 = s22 + a(m-1,l)*b(l,n-1)
406 s23 = s23 + a(m-1,l)*b(l,n)
407 s31 = s31 + a(m,l)*b(l,n-2)
408 s32 = s32 + a(m,l)*b(l,n-1)
409 s33 = s33 + a(m,l)*b(l,n)
410 enddo
411 c(m-2,n-2) = s11
412 c(m-2,n-1) = s12
413 c(m-2,n) = s13
414 c(m-1,n-2) = s21
415 c(m-1,n-1) = s22
416 c(m-1,n) = s23
417 c(m,n-2) = s31
418 c(m,n-1) = s32
419 c(m,n) = s33
420 return
421 endif
422 endif
423
424 return
425 end
Note: See TracBrowser for help on using the repository browser.