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