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