source: CIVL/examples/fortran/provesaExamples/MXM/mxm_std.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: 14.5 KB
Line 
1 subroutine mxmf2(A,N1,B,N2,C,N3)
2c
3c unrolled loop version
4c
5 real a(n1,n2),b(n2,n3),c(n1,n3)
6
7 if (n2.le.8) then
8 if (n2.eq.1) then
9 call mxf1(a,n1,b,n2,c,n3)
10 elseif (n2.eq.2) then
11 call mxf2(a,n1,b,n2,c,n3)
12 elseif (n2.eq.3) then
13 call mxf3(a,n1,b,n2,c,n3)
14 elseif (n2.eq.4) then
15 call mxf4(a,n1,b,n2,c,n3)
16 elseif (n2.eq.5) then
17 call mxf5(a,n1,b,n2,c,n3)
18 elseif (n2.eq.6) then
19 call mxf6(a,n1,b,n2,c,n3)
20 elseif (n2.eq.7) then
21 call mxf7(a,n1,b,n2,c,n3)
22 else
23 call mxf8(a,n1,b,n2,c,n3)
24 endif
25 elseif (n2.le.16) then
26 if (n2.eq.9) then
27 call mxf9(a,n1,b,n2,c,n3)
28 elseif (n2.eq.10) then
29 call mxf10(a,n1,b,n2,c,n3)
30 elseif (n2.eq.11) then
31 call mxf11(a,n1,b,n2,c,n3)
32 elseif (n2.eq.12) then
33 call mxf12(a,n1,b,n2,c,n3)
34 elseif (n2.eq.13) then
35 call mxf13(a,n1,b,n2,c,n3)
36 elseif (n2.eq.14) then
37 call mxf14(a,n1,b,n2,c,n3)
38 elseif (n2.eq.15) then
39 call mxf15(a,n1,b,n2,c,n3)
40 else
41 call mxf16(a,n1,b,n2,c,n3)
42 endif
43 elseif (n2.le.24) then
44 if (n2.eq.17) then
45 call mxf17(a,n1,b,n2,c,n3)
46 elseif (n2.eq.18) then
47 call mxf18(a,n1,b,n2,c,n3)
48 elseif (n2.eq.19) then
49 call mxf19(a,n1,b,n2,c,n3)
50 elseif (n2.eq.20) then
51 call mxf20(a,n1,b,n2,c,n3)
52 elseif (n2.eq.21) then
53 call mxf21(a,n1,b,n2,c,n3)
54 elseif (n2.eq.22) then
55 call mxf22(a,n1,b,n2,c,n3)
56 elseif (n2.eq.23) then
57 call mxf23(a,n1,b,n2,c,n3)
58 elseif (n2.eq.24) then
59 call mxf24(a,n1,b,n2,c,n3)
60 endif
61 else
62 call mxm44_0(a,n1,b,n2,c,n3)
63 endif
64
65 return
66 end
67c-----------------------------------------------------------------------
68 subroutine mxf1(a,n1,b,n2,c,n3)
69c
70 real a(n1,1),b(1,n3),c(n1,n3)
71c
72 do j=1,n3
73 do i=1,n1
74 c(i,j) = a(i,1)*b(1,j)
75 enddo
76 enddo
77 return
78 end
79c-----------------------------------------------------------------------
80 subroutine mxf2(a,n1,b,n2,c,n3)
81c
82 real a(n1,2),b(2,n3),c(n1,n3)
83c
84 do j=1,n3
85 do i=1,n1
86 c(i,j) = a(i,1)*b(1,j)
87 $ + a(i,2)*b(2,j)
88 enddo
89 enddo
90 return
91 end
92c-----------------------------------------------------------------------
93 subroutine mxf3(a,n1,b,n2,c,n3)
94c
95 real a(n1,3),b(3,n3),c(n1,n3)
96c
97 do j=1,n3
98 do i=1,n1
99 c(i,j) = a(i,1)*b(1,j)
100 $ + a(i,2)*b(2,j)
101 $ + a(i,3)*b(3,j)
102 enddo
103 enddo
104 return
105 end
106c-----------------------------------------------------------------------
107 subroutine mxf4(a,n1,b,n2,c,n3)
108c
109 real a(n1,4),b(4,n3),c(n1,n3)
110c
111 do j=1,n3
112 do i=1,n1
113 c(i,j) = a(i,1)*b(1,j)
114 $ + a(i,2)*b(2,j)
115 $ + a(i,3)*b(3,j)
116 $ + a(i,4)*b(4,j)
117 enddo
118 enddo
119 return
120 end
121c-----------------------------------------------------------------------
122c-----------------------------------------------------------------------
123c subroutine mxf5 .. mxf24 omitted for brevity
124c-----------------------------------------------------------------------
125c-----------------------------------------------------------------------
126 subroutine mxm44_0(a, m, b, k, c, n)
127c
128c matrix multiply with a 4x4 pencil
129c
130 real a(m,k), b(k,n), c(m,n)
131 real s11, s12, s13, s14, s21, s22, s23, s24
132 real s31, s32, s33, s34, s41, s42, s43, s44
133
134 mresid = iand(m,3)
135 nresid = iand(n,3)
136 m1 = m - mresid + 1
137 n1 = n - nresid + 1
138
139 do i=1,m-mresid,4
140 do j=1,n-nresid,4
141 s11 = 0.0d0
142 s21 = 0.0d0
143 s31 = 0.0d0
144 s41 = 0.0d0
145 s12 = 0.0d0
146 s22 = 0.0d0
147 s32 = 0.0d0
148 s42 = 0.0d0
149 s13 = 0.0d0
150 s23 = 0.0d0
151 s33 = 0.0d0
152 s43 = 0.0d0
153 s14 = 0.0d0
154 s24 = 0.0d0
155 s34 = 0.0d0
156 s44 = 0.0d0
157 do l=1,k
158 s11 = s11 + a(i,l)*b(l,j)
159 s12 = s12 + a(i,l)*b(l,j+1)
160 s13 = s13 + a(i,l)*b(l,j+2)
161 s14 = s14 + a(i,l)*b(l,j+3)
162
163 s21 = s21 + a(i+1,l)*b(l,j)
164 s22 = s22 + a(i+1,l)*b(l,j+1)
165 s23 = s23 + a(i+1,l)*b(l,j+2)
166 s24 = s24 + a(i+1,l)*b(l,j+3)
167
168 s31 = s31 + a(i+2,l)*b(l,j)
169 s32 = s32 + a(i+2,l)*b(l,j+1)
170 s33 = s33 + a(i+2,l)*b(l,j+2)
171 s34 = s34 + a(i+2,l)*b(l,j+3)
172
173 s41 = s41 + a(i+3,l)*b(l,j)
174 s42 = s42 + a(i+3,l)*b(l,j+1)
175 s43 = s43 + a(i+3,l)*b(l,j+2)
176 s44 = s44 + a(i+3,l)*b(l,j+3)
177 enddo
178 c(i,j) = s11
179 c(i,j+1) = s12
180 c(i,j+2) = s13
181 c(i,j+3) = s14
182
183 c(i+1,j) = s21
184 c(i+2,j) = s31
185 c(i+3,j) = s41
186
187 c(i+1,j+1) = s22
188 c(i+2,j+1) = s32
189 c(i+3,j+1) = s42
190
191 c(i+1,j+2) = s23
192 c(i+2,j+2) = s33
193 c(i+3,j+2) = s43
194
195 c(i+1,j+3) = s24
196 c(i+2,j+3) = s34
197 c(i+3,j+3) = s44
198 enddo
199* Residual when n is not multiple of 4
200 if (nresid .ne. 0) then
201 if (nresid .eq. 1) then
202 s11 = 0.0d0
203 s21 = 0.0d0
204 s31 = 0.0d0
205 s41 = 0.0d0
206 do l=1,k
207 s11 = s11 + a(i,l)*b(l,n)
208 s21 = s21 + a(i+1,l)*b(l,n)
209 s31 = s31 + a(i+2,l)*b(l,n)
210 s41 = s41 + a(i+3,l)*b(l,n)
211 enddo
212 c(i,n) = s11
213 c(i+1,n) = s21
214 c(i+2,n) = s31
215 c(i+3,n) = s41
216 elseif (nresid .eq. 2) then
217 s11 = 0.0d0
218 s21 = 0.0d0
219 s31 = 0.0d0
220 s41 = 0.0d0
221 s12 = 0.0d0
222 s22 = 0.0d0
223 s32 = 0.0d0
224 s42 = 0.0d0
225 do l=1,k
226 s11 = s11 + a(i,l)*b(l,j)
227 s12 = s12 + a(i,l)*b(l,j+1)
228
229 s21 = s21 + a(i+1,l)*b(l,j)
230 s22 = s22 + a(i+1,l)*b(l,j+1)
231
232 s31 = s31 + a(i+2,l)*b(l,j)
233 s32 = s32 + a(i+2,l)*b(l,j+1)
234
235 s41 = s41 + a(i+3,l)*b(l,j)
236 s42 = s42 + a(i+3,l)*b(l,j+1)
237 enddo
238 c(i,j) = s11
239 c(i,j+1) = s12
240
241 c(i+1,j) = s21
242 c(i+2,j) = s31
243 c(i+3,j) = s41
244
245 c(i+1,j+1) = s22
246 c(i+2,j+1) = s32
247 c(i+3,j+1) = s42
248 else
249 s11 = 0.0d0
250 s21 = 0.0d0
251 s31 = 0.0d0
252 s41 = 0.0d0
253 s12 = 0.0d0
254 s22 = 0.0d0
255 s32 = 0.0d0
256 s42 = 0.0d0
257 s13 = 0.0d0
258 s23 = 0.0d0
259 s33 = 0.0d0
260 s43 = 0.0d0
261 do l=1,k
262 s11 = s11 + a(i,l)*b(l,j)
263 s12 = s12 + a(i,l)*b(l,j+1)
264 s13 = s13 + a(i,l)*b(l,j+2)
265
266 s21 = s21 + a(i+1,l)*b(l,j)
267 s22 = s22 + a(i+1,l)*b(l,j+1)
268 s23 = s23 + a(i+1,l)*b(l,j+2)
269
270 s31 = s31 + a(i+2,l)*b(l,j)
271 s32 = s32 + a(i+2,l)*b(l,j+1)
272 s33 = s33 + a(i+2,l)*b(l,j+2)
273
274 s41 = s41 + a(i+3,l)*b(l,j)
275 s42 = s42 + a(i+3,l)*b(l,j+1)
276 s43 = s43 + a(i+3,l)*b(l,j+2)
277 enddo
278 c(i,j) = s11
279 c(i+1,j) = s21
280 c(i+2,j) = s31
281 c(i+3,j) = s41
282 c(i,j+1) = s12
283 c(i+1,j+1) = s22
284 c(i+2,j+1) = s32
285 c(i+3,j+1) = s42
286 c(i,j+2) = s13
287 c(i+1,j+2) = s23
288 c(i+2,j+2) = s33
289 c(i+3,j+2) = s43
290 endif
291 endif
292 enddo
293
294* Residual when m is not multiple of 4
295 if (mresid .eq. 0) then
296 return
297 elseif (mresid .eq. 1) then
298 do j=1,n-nresid,4
299 s11 = 0.0d0
300 s12 = 0.0d0
301 s13 = 0.0d0
302 s14 = 0.0d0
303 do l=1,k
304 s11 = s11 + a(m,l)*b(l,j)
305 s12 = s12 + a(m,l)*b(l,j+1)
306 s13 = s13 + a(m,l)*b(l,j+2)
307 s14 = s14 + a(m,l)*b(l,j+3)
308 enddo
309 c(m,j) = s11
310 c(m,j+1) = s12
311 c(m,j+2) = s13
312 c(m,j+3) = s14
313 enddo
314* mresid is 1, check nresid
315 if (nresid .eq. 0) then
316 return
317 elseif (nresid .eq. 1) then
318 s11 = 0.0d0
319 do l=1,k
320 s11 = s11 + a(m,l)*b(l,n)
321 enddo
322 c(m,n) = s11
323 return
324 elseif (nresid .eq. 2) then
325 s11 = 0.0d0
326 s12 = 0.0d0
327 do l=1,k
328 s11 = s11 + a(m,l)*b(l,n-1)
329 s12 = s12 + a(m,l)*b(l,n)
330 enddo
331 c(m,n-1) = s11
332 c(m,n) = s12
333 return
334 else
335 s11 = 0.0d0
336 s12 = 0.0d0
337 s13 = 0.0d0
338 do l=1,k
339 s11 = s11 + a(m,l)*b(l,n-2)
340 s12 = s12 + a(m,l)*b(l,n-1)
341 s13 = s13 + a(m,l)*b(l,n)
342 enddo
343 c(m,n-2) = s11
344 c(m,n-1) = s12
345 c(m,n) = s13
346 return
347 endif
348 elseif (mresid .eq. 2) then
349 do j=1,n-nresid,4
350 s11 = 0.0d0
351 s12 = 0.0d0
352 s13 = 0.0d0
353 s14 = 0.0d0
354 s21 = 0.0d0
355 s22 = 0.0d0
356 s23 = 0.0d0
357 s24 = 0.0d0
358 do l=1,k
359 s11 = s11 + a(m-1,l)*b(l,j)
360 s12 = s12 + a(m-1,l)*b(l,j+1)
361 s13 = s13 + a(m-1,l)*b(l,j+2)
362 s14 = s14 + a(m-1,l)*b(l,j+3)
363
364 s21 = s21 + a(m,l)*b(l,j)
365 s22 = s22 + a(m,l)*b(l,j+1)
366 s23 = s23 + a(m,l)*b(l,j+2)
367 s24 = s24 + a(m,l)*b(l,j+3)
368 enddo
369 c(m-1,j) = s11
370 c(m-1,j+1) = s12
371 c(m-1,j+2) = s13
372 c(m-1,j+3) = s14
373 c(m,j) = s21
374 c(m,j+1) = s22
375 c(m,j+2) = s23
376 c(m,j+3) = s24
377 enddo
378* mresid is 2, check nresid
379 if (nresid .eq. 0) then
380 return
381 elseif (nresid .eq. 1) then
382 s11 = 0.0d0
383 s21 = 0.0d0
384 do l=1,k
385 s11 = s11 + a(m-1,l)*b(l,n)
386 s21 = s21 + a(m,l)*b(l,n)
387 enddo
388 c(m-1,n) = s11
389 c(m,n) = s21
390 return
391 elseif (nresid .eq. 2) then
392 s11 = 0.0d0
393 s21 = 0.0d0
394 s12 = 0.0d0
395 s22 = 0.0d0
396 do l=1,k
397 s11 = s11 + a(m-1,l)*b(l,n-1)
398 s12 = s12 + a(m-1,l)*b(l,n)
399 s21 = s21 + a(m,l)*b(l,n-1)
400 s22 = s22 + a(m,l)*b(l,n)
401 enddo
402 c(m-1,n-1) = s11
403 c(m-1,n) = s12
404 c(m,n-1) = s21
405 c(m,n) = s22
406 return
407 else
408 s11 = 0.0d0
409 s21 = 0.0d0
410 s12 = 0.0d0
411 s22 = 0.0d0
412 s13 = 0.0d0
413 s23 = 0.0d0
414 do l=1,k
415 s11 = s11 + a(m-1,l)*b(l,n-2)
416 s12 = s12 + a(m-1,l)*b(l,n-1)
417 s13 = s13 + a(m-1,l)*b(l,n)
418 s21 = s21 + a(m,l)*b(l,n-2)
419 s22 = s22 + a(m,l)*b(l,n-1)
420 s23 = s23 + a(m,l)*b(l,n)
421 enddo
422 c(m-1,n-2) = s11
423 c(m-1,n-1) = s12
424 c(m-1,n) = s13
425 c(m,n-2) = s21
426 c(m,n-1) = s22
427 c(m,n) = s23
428 return
429 endif
430 else
431* mresid is 3
432 do j=1,n-nresid,4
433 s11 = 0.0d0
434 s21 = 0.0d0
435 s31 = 0.0d0
436
437 s12 = 0.0d0
438 s22 = 0.0d0
439 s32 = 0.0d0
440
441 s13 = 0.0d0
442 s23 = 0.0d0
443 s33 = 0.0d0
444
445 s14 = 0.0d0
446 s24 = 0.0d0
447 s34 = 0.0d0
448
449 do l=1,k
450 s11 = s11 + a(m-2,l)*b(l,j)
451 s12 = s12 + a(m-2,l)*b(l,j+1)
452 s13 = s13 + a(m-2,l)*b(l,j+2)
453 s14 = s14 + a(m-2,l)*b(l,j+3)
454
455 s21 = s21 + a(m-1,l)*b(l,j)
456 s22 = s22 + a(m-1,l)*b(l,j+1)
457 s23 = s23 + a(m-1,l)*b(l,j+2)
458 s24 = s24 + a(m-1,l)*b(l,j+3)
459
460 s31 = s31 + a(m,l)*b(l,j)
461 s32 = s32 + a(m,l)*b(l,j+1)
462 s33 = s33 + a(m,l)*b(l,j+2)
463 s34 = s34 + a(m,l)*b(l,j+3)
464 enddo
465 c(m-2,j) = s11
466 c(m-2,j+1) = s12
467 c(m-2,j+2) = s13
468 c(m-2,j+3) = s14
469
470 c(m-1,j) = s21
471 c(m-1,j+1) = s22
472 c(m-1,j+2) = s23
473 c(m-1,j+3) = s24
474
475 c(m,j) = s31
476 c(m,j+1) = s32
477 c(m,j+2) = s33
478 c(m,j+3) = s34
479 enddo
480* mresid is 3, check nresid
481 if (nresid .eq. 0) then
482 return
483 elseif (nresid .eq. 1) then
484 s11 = 0.0d0
485 s21 = 0.0d0
486 s31 = 0.0d0
487 do l=1,k
488 s11 = s11 + a(m-2,l)*b(l,n)
489 s21 = s21 + a(m-1,l)*b(l,n)
490 s31 = s31 + a(m,l)*b(l,n)
491 enddo
492 c(m-2,n) = s11
493 c(m-1,n) = s21
494 c(m,n) = s31
495 return
496 elseif (nresid .eq. 2) then
497 s11 = 0.0d0
498 s21 = 0.0d0
499 s31 = 0.0d0
500 s12 = 0.0d0
501 s22 = 0.0d0
502 s32 = 0.0d0
503 do l=1,k
504 s11 = s11 + a(m-2,l)*b(l,n-1)
505 s12 = s12 + a(m-2,l)*b(l,n)
506 s21 = s21 + a(m-1,l)*b(l,n-1)
507 s22 = s22 + a(m-1,l)*b(l,n)
508 s31 = s31 + a(m,l)*b(l,n-1)
509 s32 = s32 + a(m,l)*b(l,n)
510 enddo
511 c(m-2,n-1) = s11
512 c(m-2,n) = s12
513 c(m-1,n-1) = s21
514 c(m-1,n) = s22
515 c(m,n-1) = s31
516 c(m,n) = s32
517 return
518 else
519 s11 = 0.0d0
520 s21 = 0.0d0
521 s31 = 0.0d0
522 s12 = 0.0d0
523 s22 = 0.0d0
524 s32 = 0.0d0
525 s13 = 0.0d0
526 s23 = 0.0d0
527 s33 = 0.0d0
528 do l=1,k
529 s11 = s11 + a(m-2,l)*b(l,n-2)
530 s12 = s12 + a(m-2,l)*b(l,n-1)
531 s13 = s13 + a(m-2,l)*b(l,n)
532 s21 = s21 + a(m-1,l)*b(l,n-2)
533 s22 = s22 + a(m-1,l)*b(l,n-1)
534 s23 = s23 + a(m-1,l)*b(l,n)
535 s31 = s31 + a(m,l)*b(l,n-2)
536 s32 = s32 + a(m,l)*b(l,n-1)
537 s33 = s33 + a(m,l)*b(l,n)
538 enddo
539 c(m-2,n-2) = s11
540 c(m-2,n-1) = s12
541 c(m-2,n) = s13
542 c(m-1,n-2) = s21
543 c(m-1,n-1) = s22
544 c(m-1,n) = s23
545 c(m,n-2) = s31
546 c(m,n-1) = s32
547 c(m,n) = s33
548 return
549 endif
550 endif
551
552 return
553 end
Note: See TracBrowser for help on using the repository browser.