source: CIVL/examples/fortran/nek5000/verification/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 100755
File size: 147.2 KB
Line 
1c-----------------------------------------------------------------------
2 subroutine mxmf2(A,N1,B,N2,C,N3)
3c
4c unrolled loop version
5c
6 real a(n1,n2),b(n2,n3),c(n1,n3)
7
8 if (n2.le.8) then
9 if (n2.eq.1) then
10 call mxf1(a,n1,b,n2,c,n3)
11 elseif (n2.eq.2) then
12 call mxf2(a,n1,b,n2,c,n3)
13 elseif (n2.eq.3) then
14 call mxf3(a,n1,b,n2,c,n3)
15 elseif (n2.eq.4) then
16 call mxf4(a,n1,b,n2,c,n3)
17 elseif (n2.eq.5) then
18 call mxf5(a,n1,b,n2,c,n3)
19 elseif (n2.eq.6) then
20 call mxf6(a,n1,b,n2,c,n3)
21 elseif (n2.eq.7) then
22 call mxf7(a,n1,b,n2,c,n3)
23 else
24 call mxf8(a,n1,b,n2,c,n3)
25 endif
26 elseif (n2.le.16) then
27 if (n2.eq.9) then
28 call mxf9(a,n1,b,n2,c,n3)
29 elseif (n2.eq.10) then
30 call mxf10(a,n1,b,n2,c,n3)
31 elseif (n2.eq.11) then
32 call mxf11(a,n1,b,n2,c,n3)
33 elseif (n2.eq.12) then
34 call mxf12(a,n1,b,n2,c,n3)
35 elseif (n2.eq.13) then
36 call mxf13(a,n1,b,n2,c,n3)
37 elseif (n2.eq.14) then
38 call mxf14(a,n1,b,n2,c,n3)
39 elseif (n2.eq.15) then
40 call mxf15(a,n1,b,n2,c,n3)
41 else
42 call mxf16(a,n1,b,n2,c,n3)
43 endif
44 elseif (n2.le.24) then
45 if (n2.eq.17) then
46 call mxf17(a,n1,b,n2,c,n3)
47 elseif (n2.eq.18) then
48 call mxf18(a,n1,b,n2,c,n3)
49 elseif (n2.eq.19) then
50 call mxf19(a,n1,b,n2,c,n3)
51 elseif (n2.eq.20) then
52 call mxf20(a,n1,b,n2,c,n3)
53 elseif (n2.eq.21) then
54 call mxf21(a,n1,b,n2,c,n3)
55 elseif (n2.eq.22) then
56 call mxf22(a,n1,b,n2,c,n3)
57 elseif (n2.eq.23) then
58 call mxf23(a,n1,b,n2,c,n3)
59 elseif (n2.eq.24) then
60 call mxf24(a,n1,b,n2,c,n3)
61 endif
62 else
63 call mxm44_0(a,n1,b,n2,c,n3)
64 endif
65c
66 return
67 end
68c-----------------------------------------------------------------------
69 subroutine mxf1(a,n1,b,n2,c,n3)
70c
71 real a(n1,1),b(1,n3),c(n1,n3)
72c
73 do j=1,n3
74 do i=1,n1
75 c(i,j) = a(i,1)*b(1,j)
76 enddo
77 enddo
78 return
79 end
80c-----------------------------------------------------------------------
81 subroutine mxf2(a,n1,b,n2,c,n3)
82c
83 real a(n1,2),b(2,n3),c(n1,n3)
84c
85 do j=1,n3
86 do i=1,n1
87 c(i,j) = a(i,1)*b(1,j)
88 $ + a(i,2)*b(2,j)
89 enddo
90 enddo
91 return
92 end
93c-----------------------------------------------------------------------
94 subroutine mxf3(a,n1,b,n2,c,n3)
95c
96 real a(n1,3),b(3,n3),c(n1,n3)
97c
98 do j=1,n3
99 do i=1,n1
100 c(i,j) = a(i,1)*b(1,j)
101 $ + a(i,2)*b(2,j)
102 $ + a(i,3)*b(3,j)
103 enddo
104 enddo
105 return
106 end
107c-----------------------------------------------------------------------
108 subroutine mxf4(a,n1,b,n2,c,n3)
109c
110 real a(n1,4),b(4,n3),c(n1,n3)
111c
112 do j=1,n3
113 do i=1,n1
114 c(i,j) = a(i,1)*b(1,j)
115 $ + a(i,2)*b(2,j)
116 $ + a(i,3)*b(3,j)
117 $ + a(i,4)*b(4,j)
118 enddo
119 enddo
120 return
121 end
122c-----------------------------------------------------------------------
123 subroutine mxf5(a,n1,b,n2,c,n3)
124c
125 real a(n1,5),b(5,n3),c(n1,n3)
126c
127 do j=1,n3
128 do i=1,n1
129 c(i,j) = a(i,1)*b(1,j)
130 $ + a(i,2)*b(2,j)
131 $ + a(i,3)*b(3,j)
132 $ + a(i,4)*b(4,j)
133 $ + a(i,5)*b(5,j)
134 enddo
135 enddo
136 return
137 end
138c-----------------------------------------------------------------------
139 subroutine mxf6(a,n1,b,n2,c,n3)
140c
141 real a(n1,6),b(6,n3),c(n1,n3)
142c
143 do j=1,n3
144 do i=1,n1
145 c(i,j) = a(i,1)*b(1,j)
146 $ + a(i,2)*b(2,j)
147 $ + a(i,3)*b(3,j)
148 $ + a(i,4)*b(4,j)
149 $ + a(i,5)*b(5,j)
150 $ + a(i,6)*b(6,j)
151 enddo
152 enddo
153 return
154 end
155c-----------------------------------------------------------------------
156 subroutine mxf7(a,n1,b,n2,c,n3)
157c
158 real a(n1,7),b(7,n3),c(n1,n3)
159c
160 do j=1,n3
161 do i=1,n1
162 c(i,j) = a(i,1)*b(1,j)
163 $ + a(i,2)*b(2,j)
164 $ + a(i,3)*b(3,j)
165 $ + a(i,4)*b(4,j)
166 $ + a(i,5)*b(5,j)
167 $ + a(i,6)*b(6,j)
168 $ + a(i,7)*b(7,j)
169 enddo
170 enddo
171 return
172 end
173c-----------------------------------------------------------------------
174 subroutine mxf8(a,n1,b,n2,c,n3)
175c
176 real a(n1,8),b(8,n3),c(n1,n3)
177c
178 do j=1,n3
179 do i=1,n1
180 c(i,j) = a(i,1)*b(1,j)
181 $ + a(i,2)*b(2,j)
182 $ + a(i,3)*b(3,j)
183 $ + a(i,4)*b(4,j)
184 $ + a(i,5)*b(5,j)
185 $ + a(i,6)*b(6,j)
186 $ + a(i,7)*b(7,j)
187 $ + a(i,8)*b(8,j)
188 enddo
189 enddo
190 return
191 end
192c-----------------------------------------------------------------------
193 subroutine mxf9(a,n1,b,n2,c,n3)
194c
195 real a(n1,9),b(9,n3),c(n1,n3)
196c
197 do j=1,n3
198 do i=1,n1
199 c(i,j) = a(i,1)*b(1,j)
200 $ + a(i,2)*b(2,j)
201 $ + a(i,3)*b(3,j)
202 $ + a(i,4)*b(4,j)
203 $ + a(i,5)*b(5,j)
204 $ + a(i,6)*b(6,j)
205 $ + a(i,7)*b(7,j)
206 $ + a(i,8)*b(8,j)
207 $ + a(i,9)*b(9,j)
208 enddo
209 enddo
210 return
211 end
212c-----------------------------------------------------------------------
213 subroutine mxf10(a,n1,b,n2,c,n3)
214c
215 real a(n1,10),b(10,n3),c(n1,n3)
216c
217 do j=1,n3
218 do i=1,n1
219 c(i,j) = a(i,1)*b(1,j)
220 $ + a(i,2)*b(2,j)
221 $ + a(i,3)*b(3,j)
222 $ + a(i,4)*b(4,j)
223 $ + a(i,5)*b(5,j)
224 $ + a(i,6)*b(6,j)
225 $ + a(i,7)*b(7,j)
226 $ + a(i,8)*b(8,j)
227 $ + a(i,9)*b(9,j)
228 $ + a(i,10)*b(10,j)
229 enddo
230 enddo
231 return
232 end
233c-----------------------------------------------------------------------
234 subroutine mxf11(a,n1,b,n2,c,n3)
235c
236 real a(n1,11),b(11,n3),c(n1,n3)
237c
238 do j=1,n3
239 do i=1,n1
240 c(i,j) = a(i,1)*b(1,j)
241 $ + a(i,2)*b(2,j)
242 $ + a(i,3)*b(3,j)
243 $ + a(i,4)*b(4,j)
244 $ + a(i,5)*b(5,j)
245 $ + a(i,6)*b(6,j)
246 $ + a(i,7)*b(7,j)
247 $ + a(i,8)*b(8,j)
248 $ + a(i,9)*b(9,j)
249 $ + a(i,10)*b(10,j)
250 $ + a(i,11)*b(11,j)
251 enddo
252 enddo
253 return
254 end
255c-----------------------------------------------------------------------
256 subroutine mxf12(a,n1,b,n2,c,n3)
257c
258 real a(n1,12),b(12,n3),c(n1,n3)
259c
260 do j=1,n3
261 do i=1,n1
262 c(i,j) = a(i,1)*b(1,j)
263 $ + a(i,2)*b(2,j)
264 $ + a(i,3)*b(3,j)
265 $ + a(i,4)*b(4,j)
266 $ + a(i,5)*b(5,j)
267 $ + a(i,6)*b(6,j)
268 $ + a(i,7)*b(7,j)
269 $ + a(i,8)*b(8,j)
270 $ + a(i,9)*b(9,j)
271 $ + a(i,10)*b(10,j)
272 $ + a(i,11)*b(11,j)
273 $ + a(i,12)*b(12,j)
274 enddo
275 enddo
276 return
277 end
278c-----------------------------------------------------------------------
279 subroutine mxf13(a,n1,b,n2,c,n3)
280c
281 real a(n1,13),b(13,n3),c(n1,n3)
282c
283 do j=1,n3
284 do i=1,n1
285 c(i,j) = a(i,1)*b(1,j)
286 $ + a(i,2)*b(2,j)
287 $ + a(i,3)*b(3,j)
288 $ + a(i,4)*b(4,j)
289 $ + a(i,5)*b(5,j)
290 $ + a(i,6)*b(6,j)
291 $ + a(i,7)*b(7,j)
292 $ + a(i,8)*b(8,j)
293 $ + a(i,9)*b(9,j)
294 $ + a(i,10)*b(10,j)
295 $ + a(i,11)*b(11,j)
296 $ + a(i,12)*b(12,j)
297 $ + a(i,13)*b(13,j)
298 enddo
299 enddo
300 return
301 end
302c-----------------------------------------------------------------------
303 subroutine mxf14(a,n1,b,n2,c,n3)
304c
305 real a(n1,14),b(14,n3),c(n1,n3)
306c
307 do j=1,n3
308 do i=1,n1
309 c(i,j) = a(i,1)*b(1,j)
310 $ + a(i,2)*b(2,j)
311 $ + a(i,3)*b(3,j)
312 $ + a(i,4)*b(4,j)
313 $ + a(i,5)*b(5,j)
314 $ + a(i,6)*b(6,j)
315 $ + a(i,7)*b(7,j)
316 $ + a(i,8)*b(8,j)
317 $ + a(i,9)*b(9,j)
318 $ + a(i,10)*b(10,j)
319 $ + a(i,11)*b(11,j)
320 $ + a(i,12)*b(12,j)
321 $ + a(i,13)*b(13,j)
322 $ + a(i,14)*b(14,j)
323 enddo
324 enddo
325 return
326 end
327c-----------------------------------------------------------------------
328 subroutine mxf15(a,n1,b,n2,c,n3)
329c
330 real a(n1,15),b(15,n3),c(n1,n3)
331c
332 do j=1,n3
333 do i=1,n1
334 c(i,j) = a(i,1)*b(1,j)
335 $ + a(i,2)*b(2,j)
336 $ + a(i,3)*b(3,j)
337 $ + a(i,4)*b(4,j)
338 $ + a(i,5)*b(5,j)
339 $ + a(i,6)*b(6,j)
340 $ + a(i,7)*b(7,j)
341 $ + a(i,8)*b(8,j)
342 $ + a(i,9)*b(9,j)
343 $ + a(i,10)*b(10,j)
344 $ + a(i,11)*b(11,j)
345 $ + a(i,12)*b(12,j)
346 $ + a(i,13)*b(13,j)
347 $ + a(i,14)*b(14,j)
348 $ + a(i,15)*b(15,j)
349 enddo
350 enddo
351 return
352 end
353c-----------------------------------------------------------------------
354 subroutine mxf16(a,n1,b,n2,c,n3)
355c
356 real a(n1,16),b(16,n3),c(n1,n3)
357c
358 do j=1,n3
359 do i=1,n1
360 c(i,j) = a(i,1)*b(1,j)
361 $ + a(i,2)*b(2,j)
362 $ + a(i,3)*b(3,j)
363 $ + a(i,4)*b(4,j)
364 $ + a(i,5)*b(5,j)
365 $ + a(i,6)*b(6,j)
366 $ + a(i,7)*b(7,j)
367 $ + a(i,8)*b(8,j)
368 $ + a(i,9)*b(9,j)
369 $ + a(i,10)*b(10,j)
370 $ + a(i,11)*b(11,j)
371 $ + a(i,12)*b(12,j)
372 $ + a(i,13)*b(13,j)
373 $ + a(i,14)*b(14,j)
374 $ + a(i,15)*b(15,j)
375 $ + a(i,16)*b(16,j)
376 enddo
377 enddo
378 return
379 end
380c-----------------------------------------------------------------------
381 subroutine mxf17(a,n1,b,n2,c,n3)
382c
383 real a(n1,17),b(17,n3),c(n1,n3)
384c
385 do j=1,n3
386 do i=1,n1
387 c(i,j) = a(i,1)*b(1,j)
388 $ + a(i,2)*b(2,j)
389 $ + a(i,3)*b(3,j)
390 $ + a(i,4)*b(4,j)
391 $ + a(i,5)*b(5,j)
392 $ + a(i,6)*b(6,j)
393 $ + a(i,7)*b(7,j)
394 $ + a(i,8)*b(8,j)
395 $ + a(i,9)*b(9,j)
396 $ + a(i,10)*b(10,j)
397 $ + a(i,11)*b(11,j)
398 $ + a(i,12)*b(12,j)
399 $ + a(i,13)*b(13,j)
400 $ + a(i,14)*b(14,j)
401 $ + a(i,15)*b(15,j)
402 $ + a(i,16)*b(16,j)
403 $ + a(i,17)*b(17,j)
404 enddo
405 enddo
406 return
407 end
408c-----------------------------------------------------------------------
409 subroutine mxf18(a,n1,b,n2,c,n3)
410c
411 real a(n1,18),b(18,n3),c(n1,n3)
412c
413 do j=1,n3
414 do i=1,n1
415 c(i,j) = a(i,1)*b(1,j)
416 $ + a(i,2)*b(2,j)
417 $ + a(i,3)*b(3,j)
418 $ + a(i,4)*b(4,j)
419 $ + a(i,5)*b(5,j)
420 $ + a(i,6)*b(6,j)
421 $ + a(i,7)*b(7,j)
422 $ + a(i,8)*b(8,j)
423 $ + a(i,9)*b(9,j)
424 $ + a(i,10)*b(10,j)
425 $ + a(i,11)*b(11,j)
426 $ + a(i,12)*b(12,j)
427 $ + a(i,13)*b(13,j)
428 $ + a(i,14)*b(14,j)
429 $ + a(i,15)*b(15,j)
430 $ + a(i,16)*b(16,j)
431 $ + a(i,17)*b(17,j)
432 $ + a(i,18)*b(18,j)
433 enddo
434 enddo
435 return
436 end
437c-----------------------------------------------------------------------
438 subroutine mxf19(a,n1,b,n2,c,n3)
439c
440 real a(n1,19),b(19,n3),c(n1,n3)
441c
442 do j=1,n3
443 do i=1,n1
444 c(i,j) = a(i,1)*b(1,j)
445 $ + a(i,2)*b(2,j)
446 $ + a(i,3)*b(3,j)
447 $ + a(i,4)*b(4,j)
448 $ + a(i,5)*b(5,j)
449 $ + a(i,6)*b(6,j)
450 $ + a(i,7)*b(7,j)
451 $ + a(i,8)*b(8,j)
452 $ + a(i,9)*b(9,j)
453 $ + a(i,10)*b(10,j)
454 $ + a(i,11)*b(11,j)
455 $ + a(i,12)*b(12,j)
456 $ + a(i,13)*b(13,j)
457 $ + a(i,14)*b(14,j)
458 $ + a(i,15)*b(15,j)
459 $ + a(i,16)*b(16,j)
460 $ + a(i,17)*b(17,j)
461 $ + a(i,18)*b(18,j)
462 $ + a(i,19)*b(19,j)
463 enddo
464 enddo
465 return
466 end
467c-----------------------------------------------------------------------
468 subroutine mxf20(a,n1,b,n2,c,n3)
469c
470 real a(n1,20),b(20,n3),c(n1,n3)
471c
472 do j=1,n3
473 do i=1,n1
474 c(i,j) = a(i,1)*b(1,j)
475 $ + a(i,2)*b(2,j)
476 $ + a(i,3)*b(3,j)
477 $ + a(i,4)*b(4,j)
478 $ + a(i,5)*b(5,j)
479 $ + a(i,6)*b(6,j)
480 $ + a(i,7)*b(7,j)
481 $ + a(i,8)*b(8,j)
482 $ + a(i,9)*b(9,j)
483 $ + a(i,10)*b(10,j)
484 $ + a(i,11)*b(11,j)
485 $ + a(i,12)*b(12,j)
486 $ + a(i,13)*b(13,j)
487 $ + a(i,14)*b(14,j)
488 $ + a(i,15)*b(15,j)
489 $ + a(i,16)*b(16,j)
490 $ + a(i,17)*b(17,j)
491 $ + a(i,18)*b(18,j)
492 $ + a(i,19)*b(19,j)
493 $ + a(i,20)*b(20,j)
494 enddo
495 enddo
496 return
497 end
498c-----------------------------------------------------------------------
499 subroutine mxf21(a,n1,b,n2,c,n3)
500c
501 real a(n1,21),b(21,n3),c(n1,n3)
502c
503 do j=1,n3
504 do i=1,n1
505 c(i,j) = a(i,1)*b(1,j)
506 $ + a(i,2)*b(2,j)
507 $ + a(i,3)*b(3,j)
508 $ + a(i,4)*b(4,j)
509 $ + a(i,5)*b(5,j)
510 $ + a(i,6)*b(6,j)
511 $ + a(i,7)*b(7,j)
512 $ + a(i,8)*b(8,j)
513 $ + a(i,9)*b(9,j)
514 $ + a(i,10)*b(10,j)
515 $ + a(i,11)*b(11,j)
516 $ + a(i,12)*b(12,j)
517 $ + a(i,13)*b(13,j)
518 $ + a(i,14)*b(14,j)
519 $ + a(i,15)*b(15,j)
520 $ + a(i,16)*b(16,j)
521 $ + a(i,17)*b(17,j)
522 $ + a(i,18)*b(18,j)
523 $ + a(i,19)*b(19,j)
524 $ + a(i,20)*b(20,j)
525 $ + a(i,21)*b(21,j)
526 enddo
527 enddo
528 return
529 end
530c-----------------------------------------------------------------------
531 subroutine mxf22(a,n1,b,n2,c,n3)
532c
533 real a(n1,22),b(22,n3),c(n1,n3)
534c
535 do j=1,n3
536 do i=1,n1
537 c(i,j) = a(i,1)*b(1,j)
538 $ + a(i,2)*b(2,j)
539 $ + a(i,3)*b(3,j)
540 $ + a(i,4)*b(4,j)
541 $ + a(i,5)*b(5,j)
542 $ + a(i,6)*b(6,j)
543 $ + a(i,7)*b(7,j)
544 $ + a(i,8)*b(8,j)
545 $ + a(i,9)*b(9,j)
546 $ + a(i,10)*b(10,j)
547 $ + a(i,11)*b(11,j)
548 $ + a(i,12)*b(12,j)
549 $ + a(i,13)*b(13,j)
550 $ + a(i,14)*b(14,j)
551 $ + a(i,15)*b(15,j)
552 $ + a(i,16)*b(16,j)
553 $ + a(i,17)*b(17,j)
554 $ + a(i,18)*b(18,j)
555 $ + a(i,19)*b(19,j)
556 $ + a(i,20)*b(20,j)
557 $ + a(i,21)*b(21,j)
558 $ + a(i,22)*b(22,j)
559 enddo
560 enddo
561 return
562 end
563c-----------------------------------------------------------------------
564 subroutine mxf23(a,n1,b,n2,c,n3)
565c
566 real a(n1,23),b(23,n3),c(n1,n3)
567c
568 do j=1,n3
569 do i=1,n1
570 c(i,j) = a(i,1)*b(1,j)
571 $ + a(i,2)*b(2,j)
572 $ + a(i,3)*b(3,j)
573 $ + a(i,4)*b(4,j)
574 $ + a(i,5)*b(5,j)
575 $ + a(i,6)*b(6,j)
576 $ + a(i,7)*b(7,j)
577 $ + a(i,8)*b(8,j)
578 $ + a(i,9)*b(9,j)
579 $ + a(i,10)*b(10,j)
580 $ + a(i,11)*b(11,j)
581 $ + a(i,12)*b(12,j)
582 $ + a(i,13)*b(13,j)
583 $ + a(i,14)*b(14,j)
584 $ + a(i,15)*b(15,j)
585 $ + a(i,16)*b(16,j)
586 $ + a(i,17)*b(17,j)
587 $ + a(i,18)*b(18,j)
588 $ + a(i,19)*b(19,j)
589 $ + a(i,20)*b(20,j)
590 $ + a(i,21)*b(21,j)
591 $ + a(i,22)*b(22,j)
592 $ + a(i,23)*b(23,j)
593 enddo
594 enddo
595 return
596 end
597c-----------------------------------------------------------------------
598 subroutine mxf24(a,n1,b,n2,c,n3)
599c
600 real a(n1,24),b(24,n3),c(n1,n3)
601c
602 do j=1,n3
603 do i=1,n1
604 c(i,j) = a(i,1)*b(1,j)
605 $ + a(i,2)*b(2,j)
606 $ + a(i,3)*b(3,j)
607 $ + a(i,4)*b(4,j)
608 $ + a(i,5)*b(5,j)
609 $ + a(i,6)*b(6,j)
610 $ + a(i,7)*b(7,j)
611 $ + a(i,8)*b(8,j)
612 $ + a(i,9)*b(9,j)
613 $ + a(i,10)*b(10,j)
614 $ + a(i,11)*b(11,j)
615 $ + a(i,12)*b(12,j)
616 $ + a(i,13)*b(13,j)
617 $ + a(i,14)*b(14,j)
618 $ + a(i,15)*b(15,j)
619 $ + a(i,16)*b(16,j)
620 $ + a(i,17)*b(17,j)
621 $ + a(i,18)*b(18,j)
622 $ + a(i,19)*b(19,j)
623 $ + a(i,20)*b(20,j)
624 $ + a(i,21)*b(21,j)
625 $ + a(i,22)*b(22,j)
626 $ + a(i,23)*b(23,j)
627 $ + a(i,24)*b(24,j)
628 enddo
629 enddo
630 return
631 end
632c-----------------------------------------------------------------------
633 subroutine mxm44_0(a, m, b, k, c, n)
634c
635c matrix multiply with a 4x4 pencil
636c
637 real a(m,k), b(k,n), c(m,n)
638 real s11, s12, s13, s14, s21, s22, s23, s24
639 real s31, s32, s33, s34, s41, s42, s43, s44
640
641 mresid = iand(m,3)
642 nresid = iand(n,3)
643 m1 = m - mresid + 1
644 n1 = n - nresid + 1
645
646 do i=1,m-mresid,4
647 do j=1,n-nresid,4
648 s11 = 0.0d0
649 s21 = 0.0d0
650 s31 = 0.0d0
651 s41 = 0.0d0
652 s12 = 0.0d0
653 s22 = 0.0d0
654 s32 = 0.0d0
655 s42 = 0.0d0
656 s13 = 0.0d0
657 s23 = 0.0d0
658 s33 = 0.0d0
659 s43 = 0.0d0
660 s14 = 0.0d0
661 s24 = 0.0d0
662 s34 = 0.0d0
663 s44 = 0.0d0
664 do l=1,k
665 s11 = s11 + a(i,l)*b(l,j)
666 s12 = s12 + a(i,l)*b(l,j+1)
667 s13 = s13 + a(i,l)*b(l,j+2)
668 s14 = s14 + a(i,l)*b(l,j+3)
669
670 s21 = s21 + a(i+1,l)*b(l,j)
671 s22 = s22 + a(i+1,l)*b(l,j+1)
672 s23 = s23 + a(i+1,l)*b(l,j+2)
673 s24 = s24 + a(i+1,l)*b(l,j+3)
674
675 s31 = s31 + a(i+2,l)*b(l,j)
676 s32 = s32 + a(i+2,l)*b(l,j+1)
677 s33 = s33 + a(i+2,l)*b(l,j+2)
678 s34 = s34 + a(i+2,l)*b(l,j+3)
679
680 s41 = s41 + a(i+3,l)*b(l,j)
681 s42 = s42 + a(i+3,l)*b(l,j+1)
682 s43 = s43 + a(i+3,l)*b(l,j+2)
683 s44 = s44 + a(i+3,l)*b(l,j+3)
684 enddo
685 c(i,j) = s11
686 c(i,j+1) = s12
687 c(i,j+2) = s13
688 c(i,j+3) = s14
689
690 c(i+1,j) = s21
691 c(i+2,j) = s31
692 c(i+3,j) = s41
693
694 c(i+1,j+1) = s22
695 c(i+2,j+1) = s32
696 c(i+3,j+1) = s42
697
698 c(i+1,j+2) = s23
699 c(i+2,j+2) = s33
700 c(i+3,j+2) = s43
701
702 c(i+1,j+3) = s24
703 c(i+2,j+3) = s34
704 c(i+3,j+3) = s44
705 enddo
706* Residual when n is not multiple of 4
707 if (nresid .ne. 0) then
708 if (nresid .eq. 1) then
709 s11 = 0.0d0
710 s21 = 0.0d0
711 s31 = 0.0d0
712 s41 = 0.0d0
713 do l=1,k
714 s11 = s11 + a(i,l)*b(l,n)
715 s21 = s21 + a(i+1,l)*b(l,n)
716 s31 = s31 + a(i+2,l)*b(l,n)
717 s41 = s41 + a(i+3,l)*b(l,n)
718 enddo
719 c(i,n) = s11
720 c(i+1,n) = s21
721 c(i+2,n) = s31
722 c(i+3,n) = s41
723 elseif (nresid .eq. 2) then
724 s11 = 0.0d0
725 s21 = 0.0d0
726 s31 = 0.0d0
727 s41 = 0.0d0
728 s12 = 0.0d0
729 s22 = 0.0d0
730 s32 = 0.0d0
731 s42 = 0.0d0
732 do l=1,k
733 s11 = s11 + a(i,l)*b(l,j)
734 s12 = s12 + a(i,l)*b(l,j+1)
735
736 s21 = s21 + a(i+1,l)*b(l,j)
737 s22 = s22 + a(i+1,l)*b(l,j+1)
738
739 s31 = s31 + a(i+2,l)*b(l,j)
740 s32 = s32 + a(i+2,l)*b(l,j+1)
741
742 s41 = s41 + a(i+3,l)*b(l,j)
743 s42 = s42 + a(i+3,l)*b(l,j+1)
744 enddo
745 c(i,j) = s11
746 c(i,j+1) = s12
747
748 c(i+1,j) = s21
749 c(i+2,j) = s31
750 c(i+3,j) = s41
751
752 c(i+1,j+1) = s22
753 c(i+2,j+1) = s32
754 c(i+3,j+1) = s42
755 else
756 s11 = 0.0d0
757 s21 = 0.0d0
758 s31 = 0.0d0
759 s41 = 0.0d0
760 s12 = 0.0d0
761 s22 = 0.0d0
762 s32 = 0.0d0
763 s42 = 0.0d0
764 s13 = 0.0d0
765 s23 = 0.0d0
766 s33 = 0.0d0
767 s43 = 0.0d0
768 do l=1,k
769 s11 = s11 + a(i,l)*b(l,j)
770 s12 = s12 + a(i,l)*b(l,j+1)
771 s13 = s13 + a(i,l)*b(l,j+2)
772
773 s21 = s21 + a(i+1,l)*b(l,j)
774 s22 = s22 + a(i+1,l)*b(l,j+1)
775 s23 = s23 + a(i+1,l)*b(l,j+2)
776
777 s31 = s31 + a(i+2,l)*b(l,j)
778 s32 = s32 + a(i+2,l)*b(l,j+1)
779 s33 = s33 + a(i+2,l)*b(l,j+2)
780
781 s41 = s41 + a(i+3,l)*b(l,j)
782 s42 = s42 + a(i+3,l)*b(l,j+1)
783 s43 = s43 + a(i+3,l)*b(l,j+2)
784 enddo
785 c(i,j) = s11
786 c(i+1,j) = s21
787 c(i+2,j) = s31
788 c(i+3,j) = s41
789 c(i,j+1) = s12
790 c(i+1,j+1) = s22
791 c(i+2,j+1) = s32
792 c(i+3,j+1) = s42
793 c(i,j+2) = s13
794 c(i+1,j+2) = s23
795 c(i+2,j+2) = s33
796 c(i+3,j+2) = s43
797 endif
798 endif
799 enddo
800
801* Residual when m is not multiple of 4
802 if (mresid .eq. 0) then
803 return
804 elseif (mresid .eq. 1) then
805 do j=1,n-nresid,4
806 s11 = 0.0d0
807 s12 = 0.0d0
808 s13 = 0.0d0
809 s14 = 0.0d0
810 do l=1,k
811 s11 = s11 + a(m,l)*b(l,j)
812 s12 = s12 + a(m,l)*b(l,j+1)
813 s13 = s13 + a(m,l)*b(l,j+2)
814 s14 = s14 + a(m,l)*b(l,j+3)
815 enddo
816 c(m,j) = s11
817 c(m,j+1) = s12
818 c(m,j+2) = s13
819 c(m,j+3) = s14
820 enddo
821* mresid is 1, check nresid
822 if (nresid .eq. 0) then
823 return
824 elseif (nresid .eq. 1) then
825 s11 = 0.0d0
826 do l=1,k
827 s11 = s11 + a(m,l)*b(l,n)
828 enddo
829 c(m,n) = s11
830 return
831 elseif (nresid .eq. 2) then
832 s11 = 0.0d0
833 s12 = 0.0d0
834 do l=1,k
835 s11 = s11 + a(m,l)*b(l,n-1)
836 s12 = s12 + a(m,l)*b(l,n)
837 enddo
838 c(m,n-1) = s11
839 c(m,n) = s12
840 return
841 else
842 s11 = 0.0d0
843 s12 = 0.0d0
844 s13 = 0.0d0
845 do l=1,k
846 s11 = s11 + a(m,l)*b(l,n-2)
847 s12 = s12 + a(m,l)*b(l,n-1)
848 s13 = s13 + a(m,l)*b(l,n)
849 enddo
850 c(m,n-2) = s11
851 c(m,n-1) = s12
852 c(m,n) = s13
853 return
854 endif
855 elseif (mresid .eq. 2) then
856 do j=1,n-nresid,4
857 s11 = 0.0d0
858 s12 = 0.0d0
859 s13 = 0.0d0
860 s14 = 0.0d0
861 s21 = 0.0d0
862 s22 = 0.0d0
863 s23 = 0.0d0
864 s24 = 0.0d0
865 do l=1,k
866 s11 = s11 + a(m-1,l)*b(l,j)
867 s12 = s12 + a(m-1,l)*b(l,j+1)
868 s13 = s13 + a(m-1,l)*b(l,j+2)
869 s14 = s14 + a(m-1,l)*b(l,j+3)
870
871 s21 = s21 + a(m,l)*b(l,j)
872 s22 = s22 + a(m,l)*b(l,j+1)
873 s23 = s23 + a(m,l)*b(l,j+2)
874 s24 = s24 + a(m,l)*b(l,j+3)
875 enddo
876 c(m-1,j) = s11
877 c(m-1,j+1) = s12
878 c(m-1,j+2) = s13
879 c(m-1,j+3) = s14
880 c(m,j) = s21
881 c(m,j+1) = s22
882 c(m,j+2) = s23
883 c(m,j+3) = s24
884 enddo
885* mresid is 2, check nresid
886 if (nresid .eq. 0) then
887 return
888 elseif (nresid .eq. 1) then
889 s11 = 0.0d0
890 s21 = 0.0d0
891 do l=1,k
892 s11 = s11 + a(m-1,l)*b(l,n)
893 s21 = s21 + a(m,l)*b(l,n)
894 enddo
895 c(m-1,n) = s11
896 c(m,n) = s21
897 return
898 elseif (nresid .eq. 2) then
899 s11 = 0.0d0
900 s21 = 0.0d0
901 s12 = 0.0d0
902 s22 = 0.0d0
903 do l=1,k
904 s11 = s11 + a(m-1,l)*b(l,n-1)
905 s12 = s12 + a(m-1,l)*b(l,n)
906 s21 = s21 + a(m,l)*b(l,n-1)
907 s22 = s22 + a(m,l)*b(l,n)
908 enddo
909 c(m-1,n-1) = s11
910 c(m-1,n) = s12
911 c(m,n-1) = s21
912 c(m,n) = s22
913 return
914 else
915 s11 = 0.0d0
916 s21 = 0.0d0
917 s12 = 0.0d0
918 s22 = 0.0d0
919 s13 = 0.0d0
920 s23 = 0.0d0
921 do l=1,k
922 s11 = s11 + a(m-1,l)*b(l,n-2)
923 s12 = s12 + a(m-1,l)*b(l,n-1)
924 s13 = s13 + a(m-1,l)*b(l,n)
925 s21 = s21 + a(m,l)*b(l,n-2)
926 s22 = s22 + a(m,l)*b(l,n-1)
927 s23 = s23 + a(m,l)*b(l,n)
928 enddo
929 c(m-1,n-2) = s11
930 c(m-1,n-1) = s12
931 c(m-1,n) = s13
932 c(m,n-2) = s21
933 c(m,n-1) = s22
934 c(m,n) = s23
935 return
936 endif
937 else
938* mresid is 3
939 do j=1,n-nresid,4
940 s11 = 0.0d0
941 s21 = 0.0d0
942 s31 = 0.0d0
943
944 s12 = 0.0d0
945 s22 = 0.0d0
946 s32 = 0.0d0
947
948 s13 = 0.0d0
949 s23 = 0.0d0
950 s33 = 0.0d0
951
952 s14 = 0.0d0
953 s24 = 0.0d0
954 s34 = 0.0d0
955
956 do l=1,k
957 s11 = s11 + a(m-2,l)*b(l,j)
958 s12 = s12 + a(m-2,l)*b(l,j+1)
959 s13 = s13 + a(m-2,l)*b(l,j+2)
960 s14 = s14 + a(m-2,l)*b(l,j+3)
961
962 s21 = s21 + a(m-1,l)*b(l,j)
963 s22 = s22 + a(m-1,l)*b(l,j+1)
964 s23 = s23 + a(m-1,l)*b(l,j+2)
965 s24 = s24 + a(m-1,l)*b(l,j+3)
966
967 s31 = s31 + a(m,l)*b(l,j)
968 s32 = s32 + a(m,l)*b(l,j+1)
969 s33 = s33 + a(m,l)*b(l,j+2)
970 s34 = s34 + a(m,l)*b(l,j+3)
971 enddo
972 c(m-2,j) = s11
973 c(m-2,j+1) = s12
974 c(m-2,j+2) = s13
975 c(m-2,j+3) = s14
976
977 c(m-1,j) = s21
978 c(m-1,j+1) = s22
979 c(m-1,j+2) = s23
980 c(m-1,j+3) = s24
981
982 c(m,j) = s31
983 c(m,j+1) = s32
984 c(m,j+2) = s33
985 c(m,j+3) = s34
986 enddo
987* mresid is 3, check nresid
988 if (nresid .eq. 0) then
989 return
990 elseif (nresid .eq. 1) then
991 s11 = 0.0d0
992 s21 = 0.0d0
993 s31 = 0.0d0
994 do l=1,k
995 s11 = s11 + a(m-2,l)*b(l,n)
996 s21 = s21 + a(m-1,l)*b(l,n)
997 s31 = s31 + a(m,l)*b(l,n)
998 enddo
999 c(m-2,n) = s11
1000 c(m-1,n) = s21
1001 c(m,n) = s31
1002 return
1003 elseif (nresid .eq. 2) then
1004 s11 = 0.0d0
1005 s21 = 0.0d0
1006 s31 = 0.0d0
1007 s12 = 0.0d0
1008 s22 = 0.0d0
1009 s32 = 0.0d0
1010 do l=1,k
1011 s11 = s11 + a(m-2,l)*b(l,n-1)
1012 s12 = s12 + a(m-2,l)*b(l,n)
1013 s21 = s21 + a(m-1,l)*b(l,n-1)
1014 s22 = s22 + a(m-1,l)*b(l,n)
1015 s31 = s31 + a(m,l)*b(l,n-1)
1016 s32 = s32 + a(m,l)*b(l,n)
1017 enddo
1018 c(m-2,n-1) = s11
1019 c(m-2,n) = s12
1020 c(m-1,n-1) = s21
1021 c(m-1,n) = s22
1022 c(m,n-1) = s31
1023 c(m,n) = s32
1024 return
1025 else
1026 s11 = 0.0d0
1027 s21 = 0.0d0
1028 s31 = 0.0d0
1029 s12 = 0.0d0
1030 s22 = 0.0d0
1031 s32 = 0.0d0
1032 s13 = 0.0d0
1033 s23 = 0.0d0
1034 s33 = 0.0d0
1035 do l=1,k
1036 s11 = s11 + a(m-2,l)*b(l,n-2)
1037 s12 = s12 + a(m-2,l)*b(l,n-1)
1038 s13 = s13 + a(m-2,l)*b(l,n)
1039 s21 = s21 + a(m-1,l)*b(l,n-2)
1040 s22 = s22 + a(m-1,l)*b(l,n-1)
1041 s23 = s23 + a(m-1,l)*b(l,n)
1042 s31 = s31 + a(m,l)*b(l,n-2)
1043 s32 = s32 + a(m,l)*b(l,n-1)
1044 s33 = s33 + a(m,l)*b(l,n)
1045 enddo
1046 c(m-2,n-2) = s11
1047 c(m-2,n-1) = s12
1048 c(m-2,n) = s13
1049 c(m-1,n-2) = s21
1050 c(m-1,n-1) = s22
1051 c(m-1,n) = s23
1052 c(m,n-2) = s31
1053 c(m,n-1) = s32
1054 c(m,n) = s33
1055 return
1056 endif
1057 endif
1058
1059 return
1060 end
1061c-----------------------------------------------------------------------
1062 subroutine mxm44_2(a, m, b, k, c, n)
1063 real a(m,2), b(2,n), c(m,n)
1064
1065 nresid = iand(n,3)
1066 n1 = n - nresid + 1
1067
1068 do j=1,n-nresid,4
1069 do i=1,m
1070 c(i,j) = a(i,1)*b(1,j)
1071 $ + a(i,2)*b(2,j)
1072 c(i,j+1) = a(i,1)*b(1,j+1)
1073 $ + a(i,2)*b(2,j+1)
1074 c(i,j+2) = a(i,1)*b(1,j+2)
1075 $ + a(i,2)*b(2,j+2)
1076 c(i,j+3) = a(i,1)*b(1,j+3)
1077 $ + a(i,2)*b(2,j+3)
1078 enddo
1079 enddo
1080 if (nresid .eq. 0) then
1081 return
1082 elseif (nresid .eq. 1) then
1083 do i=1,m
1084 c(i,n) = a(i,1)*b(1,n)
1085 $ + a(i,2)*b(2,n)
1086 enddo
1087 elseif (nresid .eq. 2) then
1088 do i=1,m
1089 c(i,n-1) = a(i,1)*b(1,n-1)
1090 $ + a(i,2)*b(2,n-1)
1091 c(i,n) = a(i,1)*b(1,n)
1092 $ + a(i,2)*b(2,n)
1093 enddo
1094 else
1095 do i=1,m
1096 c(i,n-2) = a(i,1)*b(1,n-2)
1097 $ + a(i,2)*b(2,n-2)
1098 c(i,n-1) = a(i,1)*b(1,n-1)
1099 $ + a(i,2)*b(2,n-1)
1100 c(i,n) = a(i,1)*b(1,n)
1101 $ + a(i,2)*b(2,n)
1102 enddo
1103 endif
1104
1105 return
1106 end
1107c-----------------------------------------------------------------------
1108 subroutine mxm_test_all(nid,ivb)
1109c
1110c Collect matrix-matrix product statistics
1111c
1112 external mxms,mxmur2,mxmur3,mxmd,mxmfb,mxmf3,mxmu4
1113 external madd,mxm,mxm44
1114c
1115 parameter (nn=24)
1116 parameter (nt=10)
1117 character*5 c(3,nt)
1118 real s(nn,2,nt,3)
1119 real a(nn,2,nt,3)
1120
1121 call nekgsync
1122
1123 do k=1,3 ! 3 tests: N^2 x N, NxN, NxN^2
1124 call mxmtest(s(1,1, 1,k),nn,c(k, 1),mxm44 ,'mxm44',k,ivb)
1125 call mxmtest(s(1,1, 2,k),nn,c(k, 2),mxms ,' std ',k,ivb)
1126 call mxmtest(s(1,1, 3,k),nn,c(k, 3),mxmur2,'mxmu2',k,ivb)
1127 call mxmtest(s(1,1, 4,k),nn,c(k, 4),mxmur3,'mxmu3',k,ivb)
1128 call mxmtest(s(1,1, 5,k),nn,c(k, 5),mxmd ,'mxmd ',k,ivb)
1129 call mxmtest(s(1,1, 6,k),nn,c(k, 6),mxmfb ,'mxmfb',k,ivb)
1130 call mxmtest(s(1,1, 7,k),nn,c(k, 7),mxmu4 ,'mxmu4',k,ivb)
1131 call mxmtest(s(1,1, 8,k),nn,c(k, 8),mxmf3 ,'mxmf3',k,ivb)
1132 if (k.eq.2) ! Add works only for NxN case
1133 $ call mxmtest(s(1,1, 9,k),nn,c(k, 9),madd ,'madd ',k,ivb)
1134 call mxmtest(s(1,1,10,k),nn,c(k,10),mxm ,'mxm ',k,ivb)
1135 enddo
1136
1137 call nekgsync
1138 if (nid.eq.0) call mxm_analyze(s,a,nn,c,nt,ivb)
1139 call nekgsync
1140
1141 return
1142 end
1143c-----------------------------------------------------------------------
1144 subroutine initab(a,b,n)
1145 real a(1),b(1)
1146 do i=1,n-1
1147 x = i
1148 k = mod(i,19) + 2
1149 l = mod(i,17) + 5
1150 m = mod(i,31) + 3
1151 a(i) = -.25*(a(i)+a(i+1)) + (x*x + k + l)/(x*x+m)
1152 b(i) = -.25*(b(i)+b(i+1)) + (x*x + k + m)/(x*x+l)
1153 enddo
1154 a(n) = -.25*(a(n)+a(n)) + (x*x + k + l)/(x*x+m)
1155 b(n) = -.25*(b(n)+b(n)) + (x*x + k + m)/(x*x+l)
1156 return
1157 end
1158c-----------------------------------------------------------------------
1159 subroutine mxms(a,n1,b,n2,c,n3)
1160C----------------------------------------------------------------------
1161C
1162C Matrix-vector product routine.
1163C NOTE: Use assembly coded routine if available.
1164C
1165C---------------------------------------------------------------------
1166 REAL A(N1,N2),B(N2,N3),C(N1,N3)
1167C
1168 N0=N1*N3
1169 DO 10 I=1,N0
1170 C(I,1)=0.
1171 10 CONTINUE
1172 DO 100 J=1,N3
1173 DO 100 K=1,N2
1174 BB=B(K,J)
1175 DO 100 I=1,N1
1176 C(I,J)=C(I,J)+A(I,K)*BB
1177 100 CONTINUE
1178 return
1179 end
1180c-----------------------------------------------------------------------
1181 subroutine mxmu4(a,n1,b,n2,c,n3)
1182C----------------------------------------------------------------------
1183C
1184C Matrix-vector product routine.
1185C NOTE: Use assembly coded routine if available.
1186C
1187C---------------------------------------------------------------------
1188 REAL A(N1,N2),B(N2,N3),C(N1,N3)
1189C
1190 N0=N1*N3
1191 DO 10 I=1,N0
1192 C(I,1)=0.
1193 10 CONTINUE
1194 i1 = n1 - mod(n1,4) + 1
1195 DO 100 J=1,N3
1196 DO 100 K=1,N2
1197 BB=B(K,J)
1198 DO I=1,N1-3,4
1199 C(I ,J)=C(I ,J)+A(I ,K)*BB
1200 C(I+1,J)=C(I+1,J)+A(I+1,K)*BB
1201 C(I+2,J)=C(I+2,J)+A(I+2,K)*BB
1202 C(I+3,J)=C(I+3,J)+A(I+3,K)*BB
1203 enddo
1204 DO i=i1,N1
1205 C(I ,J)=C(I ,J)+A(I ,K)*BB
1206 enddo
1207 100 CONTINUE
1208 return
1209 end
1210c-----------------------------------------------------------------------
1211 subroutine madd (a,n1,b,n2,c,n3)
1212c
1213 real a(n1,n2),b(n2,n3),c(n1,n3)
1214c
1215 do j=1,n3
1216 do i=1,n1
1217 c(i,j) = a(i,j)+b(i,j)
1218 enddo
1219 enddo
1220c
1221 return
1222 end
1223c-----------------------------------------------------------------------
1224 subroutine mxmUR2(a,n1,b,n2,c,n3)
1225C----------------------------------------------------------------------
1226C
1227C Matrix-vector product routine.
1228C NOTE: Use assembly coded routine if available.
1229C
1230C---------------------------------------------------------------------
1231 REAL A(N1,N2),B(N2,N3),C(N1,N3)
1232C
1233 if (n2.le.8) then
1234 if (n2.eq.1) then
1235 call mxmur2_1(a,n1,b,n2,c,n3)
1236 elseif (n2.eq.2) then
1237 call mxmur2_2(a,n1,b,n2,c,n3)
1238 elseif (n2.eq.3) then
1239 call mxmur2_3(a,n1,b,n2,c,n3)
1240 elseif (n2.eq.4) then
1241 call mxmur2_4(a,n1,b,n2,c,n3)
1242 elseif (n2.eq.5) then
1243 call mxmur2_5(a,n1,b,n2,c,n3)
1244 elseif (n2.eq.6) then
1245 call mxmur2_6(a,n1,b,n2,c,n3)
1246 elseif (n2.eq.7) then
1247 call mxmur2_7(a,n1,b,n2,c,n3)
1248 else
1249 call mxmur2_8(a,n1,b,n2,c,n3)
1250 endif
1251 elseif (n2.le.16) then
1252 if (n2.eq.9) then
1253 call mxmur2_9(a,n1,b,n2,c,n3)
1254 elseif (n2.eq.10) then
1255 call mxmur2_10(a,n1,b,n2,c,n3)
1256 elseif (n2.eq.11) then
1257 call mxmur2_11(a,n1,b,n2,c,n3)
1258 elseif (n2.eq.12) then
1259 call mxmur2_12(a,n1,b,n2,c,n3)
1260 elseif (n2.eq.13) then
1261 call mxmur2_13(a,n1,b,n2,c,n3)
1262 elseif (n2.eq.14) then
1263 call mxmur2_14(a,n1,b,n2,c,n3)
1264 elseif (n2.eq.15) then
1265 call mxmur2_15(a,n1,b,n2,c,n3)
1266 else
1267 call mxmur2_16(a,n1,b,n2,c,n3)
1268 endif
1269 else
1270 N0=N1*N3
1271 DO 10 I=1,N0
1272 C(I,1)=0.
1273 10 CONTINUE
1274 DO 100 J=1,N3
1275 DO 100 K=1,N2
1276 BB=B(K,J)
1277 DO 100 I=1,N1
1278 C(I,J)=C(I,J)+A(I,K)*BB
1279 100 CONTINUE
1280 endif
1281 return
1282 end
1283c
1284 subroutine mxmur2_1(a,n1,b,n2,c,n3)
1285c
1286 real a(n1,1),b(1,n3),c(n1,n3)
1287c
1288 do j=1,n3
1289 do i=1,n1
1290 c(i,j) = a(i,1)*b(1,j)
1291 enddo
1292 enddo
1293 return
1294 end
1295 subroutine mxmur2_2(a,n1,b,n2,c,n3)
1296c
1297 real a(n1,2),b(2,n3),c(n1,n3)
1298c
1299 do j=1,n3
1300 do i=1,n1
1301 c(i,j) = a(i,1)*b(1,j)
1302 $ + a(i,2)*b(2,j)
1303 enddo
1304 enddo
1305 return
1306 end
1307 subroutine mxmur2_3(a,n1,b,n2,c,n3)
1308c
1309 real a(n1,3),b(3,n3),c(n1,n3)
1310c
1311 do j=1,n3
1312 do i=1,n1
1313 c(i,j) = a(i,1)*b(1,j)
1314 $ + a(i,2)*b(2,j)
1315 $ + a(i,3)*b(3,j)
1316 enddo
1317 enddo
1318 return
1319 end
1320 subroutine mxmur2_4(a,n1,b,n2,c,n3)
1321c
1322 real a(n1,4),b(4,n3),c(n1,n3)
1323c
1324 do j=1,n3
1325 do i=1,n1
1326 c(i,j) = a(i,1)*b(1,j)
1327 $ + a(i,2)*b(2,j)
1328 $ + a(i,3)*b(3,j)
1329 $ + a(i,4)*b(4,j)
1330 enddo
1331 enddo
1332 return
1333 end
1334 subroutine mxmur2_5(a,n1,b,n2,c,n3)
1335c
1336 real a(n1,5),b(5,n3),c(n1,n3)
1337c
1338 do j=1,n3
1339 do i=1,n1
1340 c(i,j) = a(i,1)*b(1,j)
1341 $ + a(i,2)*b(2,j)
1342 $ + a(i,3)*b(3,j)
1343 $ + a(i,4)*b(4,j)
1344 $ + a(i,5)*b(5,j)
1345 enddo
1346 enddo
1347 return
1348 end
1349 subroutine mxmur2_6(a,n1,b,n2,c,n3)
1350c
1351 real a(n1,6),b(6,n3),c(n1,n3)
1352c
1353 do j=1,n3
1354 do i=1,n1
1355 c(i,j) = a(i,1)*b(1,j)
1356 $ + a(i,2)*b(2,j)
1357 $ + a(i,3)*b(3,j)
1358 $ + a(i,4)*b(4,j)
1359 $ + a(i,5)*b(5,j)
1360 $ + a(i,6)*b(6,j)
1361 enddo
1362 enddo
1363 return
1364 end
1365 subroutine mxmur2_7(a,n1,b,n2,c,n3)
1366c
1367 real a(n1,7),b(7,n3),c(n1,n3)
1368c
1369 do j=1,n3
1370 do i=1,n1
1371 c(i,j) = a(i,1)*b(1,j)
1372 $ + a(i,2)*b(2,j)
1373 $ + a(i,3)*b(3,j)
1374 $ + a(i,4)*b(4,j)
1375 $ + a(i,5)*b(5,j)
1376 $ + a(i,6)*b(6,j)
1377 $ + a(i,7)*b(7,j)
1378 enddo
1379 enddo
1380 return
1381 end
1382 subroutine mxmur2_8(a,n1,b,n2,c,n3)
1383c
1384 real a(n1,8),b(8,n3),c(n1,n3)
1385c
1386 do j=1,n3
1387 do i=1,n1
1388 c(i,j) = a(i,1)*b(1,j)
1389 $ + a(i,2)*b(2,j)
1390 $ + a(i,3)*b(3,j)
1391 $ + a(i,4)*b(4,j)
1392 $ + a(i,5)*b(5,j)
1393 $ + a(i,6)*b(6,j)
1394 $ + a(i,7)*b(7,j)
1395 $ + a(i,8)*b(8,j)
1396 enddo
1397 enddo
1398 return
1399 end
1400 subroutine mxmur2_9(a,n1,b,n2,c,n3)
1401c
1402 real a(n1,9),b(9,n3),c(n1,n3)
1403c
1404 do j=1,n3
1405 do i=1,n1
1406 c(i,j) = a(i,1)*b(1,j)
1407 $ + a(i,2)*b(2,j)
1408 $ + a(i,3)*b(3,j)
1409 $ + a(i,4)*b(4,j)
1410 $ + a(i,5)*b(5,j)
1411 $ + a(i,6)*b(6,j)
1412 $ + a(i,7)*b(7,j)
1413 $ + a(i,8)*b(8,j)
1414 $ + a(i,9)*b(9,j)
1415 enddo
1416 enddo
1417 return
1418 end
1419 subroutine mxmur2_10(a,n1,b,n2,c,n3)
1420c
1421 real a(n1,10),b(10,n3),c(n1,n3)
1422c
1423 do j=1,n3
1424 do i=1,n1
1425 c(i,j) = a(i,1)*b(1,j)
1426 $ + a(i,2)*b(2,j)
1427 $ + a(i,3)*b(3,j)
1428 $ + a(i,4)*b(4,j)
1429 $ + a(i,5)*b(5,j)
1430 $ + a(i,6)*b(6,j)
1431 $ + a(i,7)*b(7,j)
1432 $ + a(i,8)*b(8,j)
1433 $ + a(i,9)*b(9,j)
1434 $ + a(i,10)*b(10,j)
1435 enddo
1436 enddo
1437 return
1438 end
1439 subroutine mxmur2_11(a,n1,b,n2,c,n3)
1440c
1441 real a(n1,11),b(11,n3),c(n1,n3)
1442c
1443 do j=1,n3
1444 do i=1,n1
1445 c(i,j) = a(i,1)*b(1,j)
1446 $ + a(i,2)*b(2,j)
1447 $ + a(i,3)*b(3,j)
1448 $ + a(i,4)*b(4,j)
1449 $ + a(i,5)*b(5,j)
1450 $ + a(i,6)*b(6,j)
1451 $ + a(i,7)*b(7,j)
1452 $ + a(i,8)*b(8,j)
1453 $ + a(i,9)*b(9,j)
1454 $ + a(i,10)*b(10,j)
1455 $ + a(i,11)*b(11,j)
1456 enddo
1457 enddo
1458 return
1459 end
1460 subroutine mxmur2_12(a,n1,b,n2,c,n3)
1461c
1462 real a(n1,12),b(12,n3),c(n1,n3)
1463c
1464 do j=1,n3
1465 do i=1,n1
1466 c(i,j) = a(i,1)*b(1,j)
1467 $ + a(i,2)*b(2,j)
1468 $ + a(i,3)*b(3,j)
1469 $ + a(i,4)*b(4,j)
1470 $ + a(i,5)*b(5,j)
1471 $ + a(i,6)*b(6,j)
1472 $ + a(i,7)*b(7,j)
1473 $ + a(i,8)*b(8,j)
1474 $ + a(i,9)*b(9,j)
1475 $ + a(i,10)*b(10,j)
1476 $ + a(i,11)*b(11,j)
1477 $ + a(i,12)*b(12,j)
1478 enddo
1479 enddo
1480 return
1481 end
1482 subroutine mxmur2_13(a,n1,b,n2,c,n3)
1483c
1484 real a(n1,13),b(13,n3),c(n1,n3)
1485c
1486 do j=1,n3
1487 do i=1,n1
1488 c(i,j) = a(i,1)*b(1,j)
1489 $ + a(i,2)*b(2,j)
1490 $ + a(i,3)*b(3,j)
1491 $ + a(i,4)*b(4,j)
1492 $ + a(i,5)*b(5,j)
1493 $ + a(i,6)*b(6,j)
1494 $ + a(i,7)*b(7,j)
1495 $ + a(i,8)*b(8,j)
1496 $ + a(i,9)*b(9,j)
1497 $ + a(i,10)*b(10,j)
1498 $ + a(i,11)*b(11,j)
1499 $ + a(i,12)*b(12,j)
1500 $ + a(i,13)*b(13,j)
1501 enddo
1502 enddo
1503 return
1504 end
1505 subroutine mxmur2_14(a,n1,b,n2,c,n3)
1506c
1507 real a(n1,14),b(14,n3),c(n1,n3)
1508c
1509 do j=1,n3
1510 do i=1,n1
1511 c(i,j) = a(i,1)*b(1,j)
1512 $ + a(i,2)*b(2,j)
1513 $ + a(i,3)*b(3,j)
1514 $ + a(i,4)*b(4,j)
1515 $ + a(i,5)*b(5,j)
1516 $ + a(i,6)*b(6,j)
1517 $ + a(i,7)*b(7,j)
1518 $ + a(i,8)*b(8,j)
1519 $ + a(i,9)*b(9,j)
1520 $ + a(i,10)*b(10,j)
1521 $ + a(i,11)*b(11,j)
1522 $ + a(i,12)*b(12,j)
1523 $ + a(i,13)*b(13,j)
1524 $ + a(i,14)*b(14,j)
1525 enddo
1526 enddo
1527 return
1528 end
1529 subroutine mxmur2_15(a,n1,b,n2,c,n3)
1530c
1531 real a(n1,15),b(15,n3),c(n1,n3)
1532c
1533 do j=1,n3
1534 do i=1,n1
1535 c(i,j) = a(i,1)*b(1,j)
1536 $ + a(i,2)*b(2,j)
1537 $ + a(i,3)*b(3,j)
1538 $ + a(i,4)*b(4,j)
1539 $ + a(i,5)*b(5,j)
1540 $ + a(i,6)*b(6,j)
1541 $ + a(i,7)*b(7,j)
1542 $ + a(i,8)*b(8,j)
1543 $ + a(i,9)*b(9,j)
1544 $ + a(i,10)*b(10,j)
1545 $ + a(i,11)*b(11,j)
1546 $ + a(i,12)*b(12,j)
1547 $ + a(i,13)*b(13,j)
1548 $ + a(i,14)*b(14,j)
1549 $ + a(i,15)*b(15,j)
1550 enddo
1551 enddo
1552 return
1553 end
1554 subroutine mxmur2_16(a,n1,b,n2,c,n3)
1555c
1556 real a(n1,16),b(16,n3),c(n1,n3)
1557c
1558 do j=1,n3
1559 do i=1,n1
1560 c(i,j) = a(i,1)*b(1,j)
1561 $ + a(i,2)*b(2,j)
1562 $ + a(i,3)*b(3,j)
1563 $ + a(i,4)*b(4,j)
1564 $ + a(i,5)*b(5,j)
1565 $ + a(i,6)*b(6,j)
1566 $ + a(i,7)*b(7,j)
1567 $ + a(i,8)*b(8,j)
1568 $ + a(i,9)*b(9,j)
1569 $ + a(i,10)*b(10,j)
1570 $ + a(i,11)*b(11,j)
1571 $ + a(i,12)*b(12,j)
1572 $ + a(i,13)*b(13,j)
1573 $ + a(i,14)*b(14,j)
1574 $ + a(i,15)*b(15,j)
1575 $ + a(i,16)*b(16,j)
1576 enddo
1577 enddo
1578 return
1579 end
1580c-----------------------------------------------------------------------
1581 subroutine mxmUR3(a,n1,b,n2,c,n3)
1582C----------------------------------------------------------------------
1583C
1584C Matrix-vector product routine.
1585C NOTE: Use assembly coded routine if available.
1586C
1587C---------------------------------------------------------------------
1588 REAL A(N1,N2),B(N2,N3),C(N1,N3)
1589C
1590 N0=N1*N3
1591 DO 10 I=1,N0
1592 C(I,1)=0.
1593 10 CONTINUE
1594 if (n3.le.8) then
1595 if (n3.eq.1) then
1596 call mxmur3_1(a,n1,b,n2,c,n3)
1597 elseif (n3.eq.2) then
1598 call mxmur3_2(a,n1,b,n2,c,n3)
1599 elseif (n3.eq.3) then
1600 call mxmur3_3(a,n1,b,n2,c,n3)
1601 elseif (n3.eq.4) then
1602 call mxmur3_4(a,n1,b,n2,c,n3)
1603 elseif (n3.eq.5) then
1604 call mxmur3_5(a,n1,b,n2,c,n3)
1605 elseif (n3.eq.6) then
1606 call mxmur3_6(a,n1,b,n2,c,n3)
1607 elseif (n3.eq.7) then
1608 call mxmur3_7(a,n1,b,n2,c,n3)
1609 else
1610 call mxmur3_8(a,n1,b,n2,c,n3)
1611 endif
1612 elseif (n3.le.16) then
1613 if (n3.eq.9) then
1614 call mxmur3_9(a,n1,b,n2,c,n3)
1615 elseif (n3.eq.10) then
1616 call mxmur3_10(a,n1,b,n2,c,n3)
1617 elseif (n3.eq.11) then
1618 call mxmur3_11(a,n1,b,n2,c,n3)
1619 elseif (n3.eq.12) then
1620 call mxmur3_12(a,n1,b,n2,c,n3)
1621 elseif (n3.eq.13) then
1622 call mxmur3_13(a,n1,b,n2,c,n3)
1623 elseif (n3.eq.14) then
1624 call mxmur3_14(a,n1,b,n2,c,n3)
1625 elseif (n3.eq.15) then
1626 call mxmur3_15(a,n1,b,n2,c,n3)
1627 else
1628 call mxmur3_16(a,n1,b,n2,c,n3)
1629 endif
1630 else
1631 DO 100 J=1,N3
1632 DO 100 K=1,N2
1633 BB=B(K,J)
1634 DO 100 I=1,N1
1635 C(I,J)=C(I,J)+A(I,K)*BB
1636 100 CONTINUE
1637 endif
1638 return
1639 end
1640c
1641 subroutine mxmur3_16(a,n1,b,n2,c,n3)
1642 real a(n1,n2),b(n2,16),c(n1,16)
1643c
1644 do k=1,n2
1645 tmp1 = b(k, 1)
1646 tmp2 = b(k, 2)
1647 tmp3 = b(k, 3)
1648 tmp4 = b(k, 4)
1649 tmp5 = b(k, 5)
1650 tmp6 = b(k, 6)
1651 tmp7 = b(k, 7)
1652 tmp8 = b(k, 8)
1653 tmp9 = b(k, 9)
1654 tmp10 = b(k,10)
1655 tmp11 = b(k,11)
1656 tmp12 = b(k,12)
1657 tmp13 = b(k,13)
1658 tmp14 = b(k,14)
1659 tmp15 = b(k,15)
1660 tmp16 = b(k,16)
1661 do i=1,n1
1662 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1663 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1664 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1665 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1666 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1667 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1668 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1669 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1670 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1671 c(i,10) = c(i,10) + a(i,k) * tmp10
1672 c(i,11) = c(i,11) + a(i,k) * tmp11
1673 c(i,12) = c(i,12) + a(i,k) * tmp12
1674 c(i,13) = c(i,13) + a(i,k) * tmp13
1675 c(i,14) = c(i,14) + a(i,k) * tmp14
1676 c(i,15) = c(i,15) + a(i,k) * tmp15
1677 c(i,16) = c(i,16) + a(i,k) * tmp16
1678 enddo
1679c
1680 enddo
1681c
1682 return
1683 end
1684 subroutine mxmur3_15(a,n1,b,n2,c,n3)
1685 real a(n1,n2),b(n2,15),c(n1,15)
1686c
1687 do k=1,n2
1688 tmp1 = b(k, 1)
1689 tmp2 = b(k, 2)
1690 tmp3 = b(k, 3)
1691 tmp4 = b(k, 4)
1692 tmp5 = b(k, 5)
1693 tmp6 = b(k, 6)
1694 tmp7 = b(k, 7)
1695 tmp8 = b(k, 8)
1696 tmp9 = b(k, 9)
1697 tmp10 = b(k,10)
1698 tmp11 = b(k,11)
1699 tmp12 = b(k,12)
1700 tmp13 = b(k,13)
1701 tmp14 = b(k,14)
1702 tmp15 = b(k,15)
1703 do i=1,n1
1704 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1705 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1706 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1707 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1708 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1709 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1710 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1711 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1712 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1713 c(i,10) = c(i,10) + a(i,k) * tmp10
1714 c(i,11) = c(i,11) + a(i,k) * tmp11
1715 c(i,12) = c(i,12) + a(i,k) * tmp12
1716 c(i,13) = c(i,13) + a(i,k) * tmp13
1717 c(i,14) = c(i,14) + a(i,k) * tmp14
1718 c(i,15) = c(i,15) + a(i,k) * tmp15
1719 enddo
1720c
1721 enddo
1722c
1723 return
1724 end
1725 subroutine mxmur3_14(a,n1,b,n2,c,n3)
1726 real a(n1,n2),b(n2,14),c(n1,14)
1727c
1728 do k=1,n2
1729 tmp1 = b(k, 1)
1730 tmp2 = b(k, 2)
1731 tmp3 = b(k, 3)
1732 tmp4 = b(k, 4)
1733 tmp5 = b(k, 5)
1734 tmp6 = b(k, 6)
1735 tmp7 = b(k, 7)
1736 tmp8 = b(k, 8)
1737 tmp9 = b(k, 9)
1738 tmp10 = b(k,10)
1739 tmp11 = b(k,11)
1740 tmp12 = b(k,12)
1741 tmp13 = b(k,13)
1742 tmp14 = b(k,14)
1743 do i=1,n1
1744 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1745 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1746 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1747 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1748 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1749 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1750 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1751 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1752 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1753 c(i,10) = c(i,10) + a(i,k) * tmp10
1754 c(i,11) = c(i,11) + a(i,k) * tmp11
1755 c(i,12) = c(i,12) + a(i,k) * tmp12
1756 c(i,13) = c(i,13) + a(i,k) * tmp13
1757 c(i,14) = c(i,14) + a(i,k) * tmp14
1758 enddo
1759c
1760 enddo
1761c
1762 return
1763 end
1764 subroutine mxmur3_13(a,n1,b,n2,c,n3)
1765 real a(n1,n2),b(n2,13),c(n1,13)
1766c
1767 do k=1,n2
1768 tmp1 = b(k, 1)
1769 tmp2 = b(k, 2)
1770 tmp3 = b(k, 3)
1771 tmp4 = b(k, 4)
1772 tmp5 = b(k, 5)
1773 tmp6 = b(k, 6)
1774 tmp7 = b(k, 7)
1775 tmp8 = b(k, 8)
1776 tmp9 = b(k, 9)
1777 tmp10 = b(k,10)
1778 tmp11 = b(k,11)
1779 tmp12 = b(k,12)
1780 tmp13 = b(k,13)
1781 do i=1,n1
1782 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1783 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1784 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1785 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1786 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1787 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1788 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1789 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1790 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1791 c(i,10) = c(i,10) + a(i,k) * tmp10
1792 c(i,11) = c(i,11) + a(i,k) * tmp11
1793 c(i,12) = c(i,12) + a(i,k) * tmp12
1794 c(i,13) = c(i,13) + a(i,k) * tmp13
1795 enddo
1796c
1797 enddo
1798c
1799 return
1800 end
1801 subroutine mxmur3_12(a,n1,b,n2,c,n3)
1802 real a(n1,n2),b(n2,12),c(n1,12)
1803c
1804 do k=1,n2
1805 tmp1 = b(k, 1)
1806 tmp2 = b(k, 2)
1807 tmp3 = b(k, 3)
1808 tmp4 = b(k, 4)
1809 tmp5 = b(k, 5)
1810 tmp6 = b(k, 6)
1811 tmp7 = b(k, 7)
1812 tmp8 = b(k, 8)
1813 tmp9 = b(k, 9)
1814 tmp10 = b(k,10)
1815 tmp11 = b(k,11)
1816 tmp12 = b(k,12)
1817 do i=1,n1
1818 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1819 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1820 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1821 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1822 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1823 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1824 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1825 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1826 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1827 c(i,10) = c(i,10) + a(i,k) * tmp10
1828 c(i,11) = c(i,11) + a(i,k) * tmp11
1829 c(i,12) = c(i,12) + a(i,k) * tmp12
1830 enddo
1831c
1832 enddo
1833c
1834 return
1835 end
1836 subroutine mxmur3_11(a,n1,b,n2,c,n3)
1837 real a(n1,n2),b(n2,11),c(n1,11)
1838c
1839 do k=1,n2
1840 tmp1 = b(k, 1)
1841 tmp2 = b(k, 2)
1842 tmp3 = b(k, 3)
1843 tmp4 = b(k, 4)
1844 tmp5 = b(k, 5)
1845 tmp6 = b(k, 6)
1846 tmp7 = b(k, 7)
1847 tmp8 = b(k, 8)
1848 tmp9 = b(k, 9)
1849 tmp10 = b(k,10)
1850 tmp11 = b(k,11)
1851 do i=1,n1
1852 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1853 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1854 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1855 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1856 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1857 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1858 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1859 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1860 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1861 c(i,10) = c(i,10) + a(i,k) * tmp10
1862 c(i,11) = c(i,11) + a(i,k) * tmp11
1863 enddo
1864c
1865 enddo
1866c
1867 return
1868 end
1869 subroutine mxmur3_10(a,n1,b,n2,c,n3)
1870 real a(n1,n2),b(n2,10),c(n1,10)
1871c
1872 do k=1,n2
1873 tmp1 = b(k, 1)
1874 tmp2 = b(k, 2)
1875 tmp3 = b(k, 3)
1876 tmp4 = b(k, 4)
1877 tmp5 = b(k, 5)
1878 tmp6 = b(k, 6)
1879 tmp7 = b(k, 7)
1880 tmp8 = b(k, 8)
1881 tmp9 = b(k, 9)
1882 tmp10 = b(k,10)
1883 do i=1,n1
1884 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1885 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1886 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1887 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1888 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1889 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1890 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1891 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1892 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1893 c(i,10) = c(i,10) + a(i,k) * tmp10
1894 enddo
1895c
1896 enddo
1897c
1898 return
1899 end
1900 subroutine mxmur3_9(a,n1,b,n2,c,n3)
1901 real a(n1,n2),b(n2,9),c(n1,9)
1902c
1903 do k=1,n2
1904 tmp1 = b(k, 1)
1905 tmp2 = b(k, 2)
1906 tmp3 = b(k, 3)
1907 tmp4 = b(k, 4)
1908 tmp5 = b(k, 5)
1909 tmp6 = b(k, 6)
1910 tmp7 = b(k, 7)
1911 tmp8 = b(k, 8)
1912 tmp9 = b(k, 9)
1913 do i=1,n1
1914 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1915 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1916 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1917 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1918 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1919 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1920 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1921 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1922 c(i, 9) = c(i, 9) + a(i,k) * tmp9
1923 enddo
1924c
1925 enddo
1926c
1927 return
1928 end
1929 subroutine mxmur3_8(a,n1,b,n2,c,n3)
1930 real a(n1,n2),b(n2,8),c(n1,8)
1931c
1932 do k=1,n2
1933 tmp1 = b(k, 1)
1934 tmp2 = b(k, 2)
1935 tmp3 = b(k, 3)
1936 tmp4 = b(k, 4)
1937 tmp5 = b(k, 5)
1938 tmp6 = b(k, 6)
1939 tmp7 = b(k, 7)
1940 tmp8 = b(k, 8)
1941 do i=1,n1
1942 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1943 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1944 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1945 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1946 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1947 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1948 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1949 c(i, 8) = c(i, 8) + a(i,k) * tmp8
1950 enddo
1951c
1952 enddo
1953c
1954 return
1955 end
1956 subroutine mxmur3_7(a,n1,b,n2,c,n3)
1957 real a(n1,n2),b(n2,7),c(n1,7)
1958c
1959 do k=1,n2
1960 tmp1 = b(k, 1)
1961 tmp2 = b(k, 2)
1962 tmp3 = b(k, 3)
1963 tmp4 = b(k, 4)
1964 tmp5 = b(k, 5)
1965 tmp6 = b(k, 6)
1966 tmp7 = b(k, 7)
1967 do i=1,n1
1968 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1969 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1970 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1971 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1972 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1973 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1974 c(i, 7) = c(i, 7) + a(i,k) * tmp7
1975 enddo
1976c
1977 enddo
1978c
1979 return
1980 end
1981 subroutine mxmur3_6(a,n1,b,n2,c,n3)
1982 real a(n1,n2),b(n2,6),c(n1,6)
1983c
1984 do k=1,n2
1985 tmp1 = b(k, 1)
1986 tmp2 = b(k, 2)
1987 tmp3 = b(k, 3)
1988 tmp4 = b(k, 4)
1989 tmp5 = b(k, 5)
1990 tmp6 = b(k, 6)
1991 do i=1,n1
1992 c(i, 1) = c(i, 1) + a(i,k) * tmp1
1993 c(i, 2) = c(i, 2) + a(i,k) * tmp2
1994 c(i, 3) = c(i, 3) + a(i,k) * tmp3
1995 c(i, 4) = c(i, 4) + a(i,k) * tmp4
1996 c(i, 5) = c(i, 5) + a(i,k) * tmp5
1997 c(i, 6) = c(i, 6) + a(i,k) * tmp6
1998 enddo
1999c
2000 enddo
2001c
2002 return
2003 end
2004 subroutine mxmur3_5(a,n1,b,n2,c,n3)
2005 real a(n1,n2),b(n2,5),c(n1,5)
2006c
2007 do k=1,n2
2008 tmp1 = b(k, 1)
2009 tmp2 = b(k, 2)
2010 tmp3 = b(k, 3)
2011 tmp4 = b(k, 4)
2012 tmp5 = b(k, 5)
2013 do i=1,n1
2014 c(i, 1) = c(i, 1) + a(i,k) * tmp1
2015 c(i, 2) = c(i, 2) + a(i,k) * tmp2
2016 c(i, 3) = c(i, 3) + a(i,k) * tmp3
2017 c(i, 4) = c(i, 4) + a(i,k) * tmp4
2018 c(i, 5) = c(i, 5) + a(i,k) * tmp5
2019 enddo
2020c
2021 enddo
2022c
2023 return
2024 end
2025 subroutine mxmur3_4(a,n1,b,n2,c,n3)
2026 real a(n1,n2),b(n2,4),c(n1,4)
2027c
2028 do k=1,n2
2029 tmp1 = b(k, 1)
2030 tmp2 = b(k, 2)
2031 tmp3 = b(k, 3)
2032 tmp4 = b(k, 4)
2033 do i=1,n1
2034 c(i, 1) = c(i, 1) + a(i,k) * tmp1
2035 c(i, 2) = c(i, 2) + a(i,k) * tmp2
2036 c(i, 3) = c(i, 3) + a(i,k) * tmp3
2037 c(i, 4) = c(i, 4) + a(i,k) * tmp4
2038 enddo
2039c
2040 enddo
2041c
2042 return
2043 end
2044 subroutine mxmur3_3(a,n1,b,n2,c,n3)
2045 real a(n1,n2),b(n2,3),c(n1,3)
2046c
2047 do k=1,n2
2048 tmp1 = b(k, 1)
2049 tmp2 = b(k, 2)
2050 tmp3 = b(k, 3)
2051 do i=1,n1
2052 c(i, 1) = c(i, 1) + a(i,k) * tmp1
2053 c(i, 2) = c(i, 2) + a(i,k) * tmp2
2054 c(i, 3) = c(i, 3) + a(i,k) * tmp3
2055 enddo
2056c
2057 enddo
2058c
2059 return
2060 end
2061 subroutine mxmur3_2(a,n1,b,n2,c,n3)
2062 real a(n1,n2),b(n2,2),c(n1,2)
2063c
2064 do k=1,n2
2065 tmp1 = b(k, 1)
2066 tmp2 = b(k, 2)
2067 do i=1,n1
2068 c(i, 1) = c(i, 1) + a(i,k) * tmp1
2069 c(i, 2) = c(i, 2) + a(i,k) * tmp2
2070 enddo
2071c
2072 enddo
2073c
2074 return
2075 end
2076 subroutine mxmur3_1(a,n1,b,n2,c,n3)
2077 real a(n1,n2),b(n2,1),c(n1,1)
2078c
2079 do k=1,n2
2080 tmp1 = b(k, 1)
2081 do i=1,n1
2082 c(i, 1) = c(i, 1) + a(i,k) * tmp1
2083 enddo
2084 enddo
2085c
2086 return
2087 end
2088C----------------------------------------------------------------------
2089 subroutine mxmd(a,n1,b,n2,c,n3)
2090C
2091C Matrix-vector product routine.
2092C NOTE: Use assembly coded routine if available.
2093C
2094C---------------------------------------------------------------------
2095 REAL A(N1,N2),B(N2,N3),C(N1,N3)
2096 REAL ONE,ZERO,EPS
2097C
2098C
2099C
2100 one=1.0
2101 zero=0.0
2102 call dgemm( 'N','N',n1,n3,n2,ONE,A,N1,B,N2,ZERO,C,N1)
2103 return
2104 end
2105c-----------------------------------------------------------------------
2106 subroutine mxmfb(a,n1,b,n2,c,n3)
2107C-----------------------------------------------------------------------
2108C
2109C Matrix-vector product routine.
2110C NOTE: Use assembly coded routine if available.
2111C
2112C----------------------------------------------------------------------
2113 REAL A(N1,N2),B(N2,N3),C(N1,N3)
2114C
2115 integer wdsize
2116 save wdsize
2117 data wdsize/0/
2118c
2119c First call: determine word size for dgemm/sgemm discrimination, below.
2120c
2121 if (wdsize.eq.0) then
2122 one = 1.0
2123 eps = 1.e-12
2124 wdsize = 8
2125 if (one+eps.eq.1.0) wdsize = 4
2126 endif
2127c
2128 if (n2.le.8) then
2129 if (n2.eq.1) then
2130 call mxmfb_1(a,n1,b,n2,c,n3)
2131 elseif (n2.eq.2) then
2132 call mxmfb_2(a,n1,b,n2,c,n3)
2133 elseif (n2.eq.3) then
2134 call mxmfb_3(a,n1,b,n2,c,n3)
2135 elseif (n2.eq.4) then
2136 call mxmfb_4(a,n1,b,n2,c,n3)
2137 elseif (n2.eq.5) then
2138 call mxmfb_5(a,n1,b,n2,c,n3)
2139 elseif (n2.eq.6) then
2140 call mxmfb_6(a,n1,b,n2,c,n3)
2141 elseif (n2.eq.7) then
2142 call mxmfb_7(a,n1,b,n2,c,n3)
2143 else
2144 call mxmfb_8(a,n1,b,n2,c,n3)
2145 endif
2146 elseif (n2.le.16) then
2147 if (n2.eq.9) then
2148 call mxmfb_9(a,n1,b,n2,c,n3)
2149 elseif (n2.eq.10) then
2150 call mxmfb_10(a,n1,b,n2,c,n3)
2151 elseif (n2.eq.11) then
2152 call mxmfb_11(a,n1,b,n2,c,n3)
2153 elseif (n2.eq.12) then
2154 call mxmfb_12(a,n1,b,n2,c,n3)
2155 elseif (n2.eq.13) then
2156 call mxmfb_13(a,n1,b,n2,c,n3)
2157 elseif (n2.eq.14) then
2158 call mxmfb_14(a,n1,b,n2,c,n3)
2159 elseif (n2.eq.15) then
2160 call mxmfb_15(a,n1,b,n2,c,n3)
2161 else
2162 call mxmfb_16(a,n1,b,n2,c,n3)
2163 endif
2164 elseif (n2.le.24) then
2165 if (n2.eq.17) then
2166 call mxmfb_17(a,n1,b,n2,c,n3)
2167 elseif (n2.eq.18) then
2168 call mxmfb_18(a,n1,b,n2,c,n3)
2169 elseif (n2.eq.19) then
2170 call mxmfb_19(a,n1,b,n2,c,n3)
2171 elseif (n2.eq.20) then
2172 call mxmfb_20(a,n1,b,n2,c,n3)
2173 elseif (n2.eq.21) then
2174 call mxmfb_21(a,n1,b,n2,c,n3)
2175 elseif (n2.eq.22) then
2176 call mxmfb_22(a,n1,b,n2,c,n3)
2177 elseif (n2.eq.23) then
2178 call mxmfb_23(a,n1,b,n2,c,n3)
2179 elseif (n2.eq.24) then
2180 call mxmfb_24(a,n1,b,n2,c,n3)
2181 endif
2182 else
2183c
2184 one=1.0
2185 zero=0.0
2186 call dgemm( 'N','N',n1,n3,n2,ONE,A,N1,B,N2,ZERO,C,N1)
2187
2188 endif
2189 return
2190 end
2191c-----------------------------------------------------------------------
2192 subroutine mxmfb_1(a,n1,b,n2,c,n3)
2193c
2194 real a(n1,1),b(1,n3),c(n1,n3)
2195c
2196 do j=1,n3
2197 do i=1,n1
2198 c(i,j) = a(i,1)*b(1,j)
2199 enddo
2200 enddo
2201 return
2202 end
2203c-----------------------------------------------------------------------
2204 subroutine mxmfb_2(a,n1,b,n2,c,n3)
2205c
2206 real a(n1,2),b(2,n3),c(n1,n3)
2207c
2208 do j=1,n3
2209 do i=1,n1
2210 c(i,j) = a(i,1)*b(1,j)
2211 $ + a(i,2)*b(2,j)
2212 enddo
2213 enddo
2214 return
2215 end
2216c-----------------------------------------------------------------------
2217 subroutine mxmfb_3(a,n1,b,n2,c,n3)
2218c
2219 real a(n1,3),b(3,n3),c(n1,n3)
2220c
2221 do j=1,n3
2222 do i=1,n1
2223 c(i,j) = a(i,1)*b(1,j)
2224 $ + a(i,2)*b(2,j)
2225 $ + a(i,3)*b(3,j)
2226 enddo
2227 enddo
2228 return
2229 end
2230c-----------------------------------------------------------------------
2231 subroutine mxmfb_4(a,n1,b,n2,c,n3)
2232c
2233 real a(n1,4),b(4,n3),c(n1,n3)
2234c
2235 do j=1,n3
2236 do i=1,n1
2237 c(i,j) = a(i,1)*b(1,j)
2238 $ + a(i,2)*b(2,j)
2239 $ + a(i,3)*b(3,j)
2240 $ + a(i,4)*b(4,j)
2241 enddo
2242 enddo
2243 return
2244 end
2245c-----------------------------------------------------------------------
2246 subroutine mxmfb_5(a,n1,b,n2,c,n3)
2247c
2248 real a(n1,5),b(5,n3),c(n1,n3)
2249c
2250 do j=1,n3
2251 do i=1,n1
2252 c(i,j) = a(i,1)*b(1,j)
2253 $ + a(i,2)*b(2,j)
2254 $ + a(i,3)*b(3,j)
2255 $ + a(i,4)*b(4,j)
2256 $ + a(i,5)*b(5,j)
2257 enddo
2258 enddo
2259 return
2260 end
2261c-----------------------------------------------------------------------
2262 subroutine mxmfb_6(a,n1,b,n2,c,n3)
2263c
2264 real a(n1,6),b(6,n3),c(n1,n3)
2265c
2266 do j=1,n3
2267 do i=1,n1
2268 c(i,j) = a(i,1)*b(1,j)
2269 $ + a(i,2)*b(2,j)
2270 $ + a(i,3)*b(3,j)
2271 $ + a(i,4)*b(4,j)
2272 $ + a(i,5)*b(5,j)
2273 $ + a(i,6)*b(6,j)
2274 enddo
2275 enddo
2276 return
2277 end
2278c-----------------------------------------------------------------------
2279 subroutine mxmfb_7(a,n1,b,n2,c,n3)
2280c
2281 real a(n1,7),b(7,n3),c(n1,n3)
2282c
2283 do j=1,n3
2284 do i=1,n1
2285 c(i,j) = a(i,1)*b(1,j)
2286 $ + a(i,2)*b(2,j)
2287 $ + a(i,3)*b(3,j)
2288 $ + a(i,4)*b(4,j)
2289 $ + a(i,5)*b(5,j)
2290 $ + a(i,6)*b(6,j)
2291 $ + a(i,7)*b(7,j)
2292 enddo
2293 enddo
2294 return
2295 end
2296c-----------------------------------------------------------------------
2297 subroutine mxmfb_8(a,n1,b,n2,c,n3)
2298c
2299 real a(n1,8),b(8,n3),c(n1,n3)
2300c
2301 do j=1,n3
2302 do i=1,n1
2303 c(i,j) = a(i,1)*b(1,j)
2304 $ + a(i,2)*b(2,j)
2305 $ + a(i,3)*b(3,j)
2306 $ + a(i,4)*b(4,j)
2307 $ + a(i,5)*b(5,j)
2308 $ + a(i,6)*b(6,j)
2309 $ + a(i,7)*b(7,j)
2310 $ + a(i,8)*b(8,j)
2311 enddo
2312 enddo
2313 return
2314 end
2315c-----------------------------------------------------------------------
2316 subroutine mxmfb_9(a,n1,b,n2,c,n3)
2317c
2318 real a(n1,9),b(9,n3),c(n1,n3)
2319c
2320 do j=1,n3
2321 do i=1,n1
2322 c(i,j) = a(i,1)*b(1,j)
2323 $ + a(i,2)*b(2,j)
2324 $ + a(i,3)*b(3,j)
2325 $ + a(i,4)*b(4,j)
2326 $ + a(i,5)*b(5,j)
2327 $ + a(i,6)*b(6,j)
2328 $ + a(i,7)*b(7,j)
2329 $ + a(i,8)*b(8,j)
2330 $ + a(i,9)*b(9,j)
2331 enddo
2332 enddo
2333 return
2334 end
2335c-----------------------------------------------------------------------
2336 subroutine mxmfb_10(a,n1,b,n2,c,n3)
2337c
2338 real a(n1,10),b(10,n3),c(n1,n3)
2339c
2340 do j=1,n3
2341 do i=1,n1
2342 c(i,j) = a(i,1)*b(1,j)
2343 $ + a(i,2)*b(2,j)
2344 $ + a(i,3)*b(3,j)
2345 $ + a(i,4)*b(4,j)
2346 $ + a(i,5)*b(5,j)
2347 $ + a(i,6)*b(6,j)
2348 $ + a(i,7)*b(7,j)
2349 $ + a(i,8)*b(8,j)
2350 $ + a(i,9)*b(9,j)
2351 $ + a(i,10)*b(10,j)
2352 enddo
2353 enddo
2354 return
2355 end
2356c-----------------------------------------------------------------------
2357 subroutine mxmfb_11(a,n1,b,n2,c,n3)
2358c
2359 real a(n1,11),b(11,n3),c(n1,n3)
2360c
2361 do j=1,n3
2362 do i=1,n1
2363 c(i,j) = a(i,1)*b(1,j)
2364 $ + a(i,2)*b(2,j)
2365 $ + a(i,3)*b(3,j)
2366 $ + a(i,4)*b(4,j)
2367 $ + a(i,5)*b(5,j)
2368 $ + a(i,6)*b(6,j)
2369 $ + a(i,7)*b(7,j)
2370 $ + a(i,8)*b(8,j)
2371 $ + a(i,9)*b(9,j)
2372 $ + a(i,10)*b(10,j)
2373 $ + a(i,11)*b(11,j)
2374 enddo
2375 enddo
2376 return
2377 end
2378c-----------------------------------------------------------------------
2379 subroutine mxmfb_12(a,n1,b,n2,c,n3)
2380c
2381 real a(n1,12),b(12,n3),c(n1,n3)
2382c
2383 do j=1,n3
2384 do i=1,n1
2385 c(i,j) = a(i,1)*b(1,j)
2386 $ + a(i,2)*b(2,j)
2387 $ + a(i,3)*b(3,j)
2388 $ + a(i,4)*b(4,j)
2389 $ + a(i,5)*b(5,j)
2390 $ + a(i,6)*b(6,j)
2391 $ + a(i,7)*b(7,j)
2392 $ + a(i,8)*b(8,j)
2393 $ + a(i,9)*b(9,j)
2394 $ + a(i,10)*b(10,j)
2395 $ + a(i,11)*b(11,j)
2396 $ + a(i,12)*b(12,j)
2397 enddo
2398 enddo
2399 return
2400 end
2401c-----------------------------------------------------------------------
2402 subroutine mxmfb_13(a,n1,b,n2,c,n3)
2403c
2404 real a(n1,13),b(13,n3),c(n1,n3)
2405c
2406 do j=1,n3
2407 do i=1,n1
2408 c(i,j) = a(i,1)*b(1,j)
2409 $ + a(i,2)*b(2,j)
2410 $ + a(i,3)*b(3,j)
2411 $ + a(i,4)*b(4,j)
2412 $ + a(i,5)*b(5,j)
2413 $ + a(i,6)*b(6,j)
2414 $ + a(i,7)*b(7,j)
2415 $ + a(i,8)*b(8,j)
2416 $ + a(i,9)*b(9,j)
2417 $ + a(i,10)*b(10,j)
2418 $ + a(i,11)*b(11,j)
2419 $ + a(i,12)*b(12,j)
2420 $ + a(i,13)*b(13,j)
2421 enddo
2422 enddo
2423 return
2424 end
2425c-----------------------------------------------------------------------
2426 subroutine mxmfb_14(a,n1,b,n2,c,n3)
2427c
2428 real a(n1,14),b(14,n3),c(n1,n3)
2429c
2430 do j=1,n3
2431 do i=1,n1
2432 c(i,j) = a(i,1)*b(1,j)
2433 $ + a(i,2)*b(2,j)
2434 $ + a(i,3)*b(3,j)
2435 $ + a(i,4)*b(4,j)
2436 $ + a(i,5)*b(5,j)
2437 $ + a(i,6)*b(6,j)
2438 $ + a(i,7)*b(7,j)
2439 $ + a(i,8)*b(8,j)
2440 $ + a(i,9)*b(9,j)
2441 $ + a(i,10)*b(10,j)
2442 $ + a(i,11)*b(11,j)
2443 $ + a(i,12)*b(12,j)
2444 $ + a(i,13)*b(13,j)
2445 $ + a(i,14)*b(14,j)
2446 enddo
2447 enddo
2448 return
2449 end
2450c-----------------------------------------------------------------------
2451 subroutine mxmfb_15(a,n1,b,n2,c,n3)
2452c
2453 real a(n1,15),b(15,n3),c(n1,n3)
2454c
2455 do j=1,n3
2456 do i=1,n1
2457 c(i,j) = a(i,1)*b(1,j)
2458 $ + a(i,2)*b(2,j)
2459 $ + a(i,3)*b(3,j)
2460 $ + a(i,4)*b(4,j)
2461 $ + a(i,5)*b(5,j)
2462 $ + a(i,6)*b(6,j)
2463 $ + a(i,7)*b(7,j)
2464 $ + a(i,8)*b(8,j)
2465 $ + a(i,9)*b(9,j)
2466 $ + a(i,10)*b(10,j)
2467 $ + a(i,11)*b(11,j)
2468 $ + a(i,12)*b(12,j)
2469 $ + a(i,13)*b(13,j)
2470 $ + a(i,14)*b(14,j)
2471 $ + a(i,15)*b(15,j)
2472 enddo
2473 enddo
2474 return
2475 end
2476c-----------------------------------------------------------------------
2477 subroutine mxmfb_16(a,n1,b,n2,c,n3)
2478c
2479 real a(n1,16),b(16,n3),c(n1,n3)
2480c
2481 do j=1,n3
2482 do i=1,n1
2483 c(i,j) = a(i,1)*b(1,j)
2484 $ + a(i,2)*b(2,j)
2485 $ + a(i,3)*b(3,j)
2486 $ + a(i,4)*b(4,j)
2487 $ + a(i,5)*b(5,j)
2488 $ + a(i,6)*b(6,j)
2489 $ + a(i,7)*b(7,j)
2490 $ + a(i,8)*b(8,j)
2491 $ + a(i,9)*b(9,j)
2492 $ + a(i,10)*b(10,j)
2493 $ + a(i,11)*b(11,j)
2494 $ + a(i,12)*b(12,j)
2495 $ + a(i,13)*b(13,j)
2496 $ + a(i,14)*b(14,j)
2497 $ + a(i,15)*b(15,j)
2498 $ + a(i,16)*b(16,j)
2499 enddo
2500 enddo
2501 return
2502 end
2503c-----------------------------------------------------------------------
2504 subroutine mxmfb_17(a,n1,b,n2,c,n3)
2505c
2506 real a(n1,17),b(17,n3),c(n1,n3)
2507c
2508 do j=1,n3
2509 do i=1,n1
2510 c(i,j) = a(i,1)*b(1,j)
2511 $ + a(i,2)*b(2,j)
2512 $ + a(i,3)*b(3,j)
2513 $ + a(i,4)*b(4,j)
2514 $ + a(i,5)*b(5,j)
2515 $ + a(i,6)*b(6,j)
2516 $ + a(i,7)*b(7,j)
2517 $ + a(i,8)*b(8,j)
2518 $ + a(i,9)*b(9,j)
2519 $ + a(i,10)*b(10,j)
2520 $ + a(i,11)*b(11,j)
2521 $ + a(i,12)*b(12,j)
2522 $ + a(i,13)*b(13,j)
2523 $ + a(i,14)*b(14,j)
2524 $ + a(i,15)*b(15,j)
2525 $ + a(i,16)*b(16,j)
2526 $ + a(i,17)*b(17,j)
2527 enddo
2528 enddo
2529 return
2530 end
2531c-----------------------------------------------------------------------
2532 subroutine mxmfb_18(a,n1,b,n2,c,n3)
2533c
2534 real a(n1,18),b(18,n3),c(n1,n3)
2535c
2536 do j=1,n3
2537 do i=1,n1
2538 c(i,j) = a(i,1)*b(1,j)
2539 $ + a(i,2)*b(2,j)
2540 $ + a(i,3)*b(3,j)
2541 $ + a(i,4)*b(4,j)
2542 $ + a(i,5)*b(5,j)
2543 $ + a(i,6)*b(6,j)
2544 $ + a(i,7)*b(7,j)
2545 $ + a(i,8)*b(8,j)
2546 $ + a(i,9)*b(9,j)
2547 $ + a(i,10)*b(10,j)
2548 $ + a(i,11)*b(11,j)
2549 $ + a(i,12)*b(12,j)
2550 $ + a(i,13)*b(13,j)
2551 $ + a(i,14)*b(14,j)
2552 $ + a(i,15)*b(15,j)
2553 $ + a(i,16)*b(16,j)
2554 $ + a(i,17)*b(17,j)
2555 $ + a(i,18)*b(18,j)
2556 enddo
2557 enddo
2558 return
2559 end
2560c-----------------------------------------------------------------------
2561 subroutine mxmfb_19(a,n1,b,n2,c,n3)
2562c
2563 real a(n1,19),b(19,n3),c(n1,n3)
2564c
2565 do j=1,n3
2566 do i=1,n1
2567 c(i,j) = a(i,1)*b(1,j)
2568 $ + a(i,2)*b(2,j)
2569 $ + a(i,3)*b(3,j)
2570 $ + a(i,4)*b(4,j)
2571 $ + a(i,5)*b(5,j)
2572 $ + a(i,6)*b(6,j)
2573 $ + a(i,7)*b(7,j)
2574 $ + a(i,8)*b(8,j)
2575 $ + a(i,9)*b(9,j)
2576 $ + a(i,10)*b(10,j)
2577 $ + a(i,11)*b(11,j)
2578 $ + a(i,12)*b(12,j)
2579 $ + a(i,13)*b(13,j)
2580 $ + a(i,14)*b(14,j)
2581 $ + a(i,15)*b(15,j)
2582 $ + a(i,16)*b(16,j)
2583 $ + a(i,17)*b(17,j)
2584 $ + a(i,18)*b(18,j)
2585 $ + a(i,19)*b(19,j)
2586 enddo
2587 enddo
2588 return
2589 end
2590c-----------------------------------------------------------------------
2591 subroutine mxmfb_20(a,n1,b,n2,c,n3)
2592c
2593 real a(n1,20),b(20,n3),c(n1,n3)
2594c
2595 do j=1,n3
2596 do i=1,n1
2597 c(i,j) = a(i,1)*b(1,j)
2598 $ + a(i,2)*b(2,j)
2599 $ + a(i,3)*b(3,j)
2600 $ + a(i,4)*b(4,j)
2601 $ + a(i,5)*b(5,j)
2602 $ + a(i,6)*b(6,j)
2603 $ + a(i,7)*b(7,j)
2604 $ + a(i,8)*b(8,j)
2605 $ + a(i,9)*b(9,j)
2606 $ + a(i,10)*b(10,j)
2607 $ + a(i,11)*b(11,j)
2608 $ + a(i,12)*b(12,j)
2609 $ + a(i,13)*b(13,j)
2610 $ + a(i,14)*b(14,j)
2611 $ + a(i,15)*b(15,j)
2612 $ + a(i,16)*b(16,j)
2613 $ + a(i,17)*b(17,j)
2614 $ + a(i,18)*b(18,j)
2615 $ + a(i,19)*b(19,j)
2616 $ + a(i,20)*b(20,j)
2617 enddo
2618 enddo
2619 return
2620 end
2621c-----------------------------------------------------------------------
2622 subroutine mxmfb_21(a,n1,b,n2,c,n3)
2623c
2624 real a(n1,21),b(21,n3),c(n1,n3)
2625c
2626 do j=1,n3
2627 do i=1,n1
2628 c(i,j) = a(i,1)*b(1,j)
2629 $ + a(i,2)*b(2,j)
2630 $ + a(i,3)*b(3,j)
2631 $ + a(i,4)*b(4,j)
2632 $ + a(i,5)*b(5,j)
2633 $ + a(i,6)*b(6,j)
2634 $ + a(i,7)*b(7,j)
2635 $ + a(i,8)*b(8,j)
2636 $ + a(i,9)*b(9,j)
2637 $ + a(i,10)*b(10,j)
2638 $ + a(i,11)*b(11,j)
2639 $ + a(i,12)*b(12,j)
2640 $ + a(i,13)*b(13,j)
2641 $ + a(i,14)*b(14,j)
2642 $ + a(i,15)*b(15,j)
2643 $ + a(i,16)*b(16,j)
2644 $ + a(i,17)*b(17,j)
2645 $ + a(i,18)*b(18,j)
2646 $ + a(i,19)*b(19,j)
2647 $ + a(i,20)*b(20,j)
2648 $ + a(i,21)*b(21,j)
2649 enddo
2650 enddo
2651 return
2652 end
2653c-----------------------------------------------------------------------
2654 subroutine mxmfb_22(a,n1,b,n2,c,n3)
2655c
2656 real a(n1,22),b(22,n3),c(n1,n3)
2657c
2658 do j=1,n3
2659 do i=1,n1
2660 c(i,j) = a(i,1)*b(1,j)
2661 $ + a(i,2)*b(2,j)
2662 $ + a(i,3)*b(3,j)
2663 $ + a(i,4)*b(4,j)
2664 $ + a(i,5)*b(5,j)
2665 $ + a(i,6)*b(6,j)
2666 $ + a(i,7)*b(7,j)
2667 $ + a(i,8)*b(8,j)
2668 $ + a(i,9)*b(9,j)
2669 $ + a(i,10)*b(10,j)
2670 $ + a(i,11)*b(11,j)
2671 $ + a(i,12)*b(12,j)
2672 $ + a(i,13)*b(13,j)
2673 $ + a(i,14)*b(14,j)
2674 $ + a(i,15)*b(15,j)
2675 $ + a(i,16)*b(16,j)
2676 $ + a(i,17)*b(17,j)
2677 $ + a(i,18)*b(18,j)
2678 $ + a(i,19)*b(19,j)
2679 $ + a(i,20)*b(20,j)
2680 $ + a(i,21)*b(21,j)
2681 $ + a(i,22)*b(22,j)
2682 enddo
2683 enddo
2684 return
2685 end
2686c-----------------------------------------------------------------------
2687 subroutine mxmfb_23(a,n1,b,n2,c,n3)
2688c
2689 real a(n1,23),b(23,n3),c(n1,n3)
2690c
2691 do j=1,n3
2692 do i=1,n1
2693 c(i,j) = a(i,1)*b(1,j)
2694 $ + a(i,2)*b(2,j)
2695 $ + a(i,3)*b(3,j)
2696 $ + a(i,4)*b(4,j)
2697 $ + a(i,5)*b(5,j)
2698 $ + a(i,6)*b(6,j)
2699 $ + a(i,7)*b(7,j)
2700 $ + a(i,8)*b(8,j)
2701 $ + a(i,9)*b(9,j)
2702 $ + a(i,10)*b(10,j)
2703 $ + a(i,11)*b(11,j)
2704 $ + a(i,12)*b(12,j)
2705 $ + a(i,13)*b(13,j)
2706 $ + a(i,14)*b(14,j)
2707 $ + a(i,15)*b(15,j)
2708 $ + a(i,16)*b(16,j)
2709 $ + a(i,17)*b(17,j)
2710 $ + a(i,18)*b(18,j)
2711 $ + a(i,19)*b(19,j)
2712 $ + a(i,20)*b(20,j)
2713 $ + a(i,21)*b(21,j)
2714 $ + a(i,22)*b(22,j)
2715 $ + a(i,23)*b(23,j)
2716 enddo
2717 enddo
2718 return
2719 end
2720c-----------------------------------------------------------------------
2721 subroutine mxmfb_24(a,n1,b,n2,c,n3)
2722c
2723 real a(n1,24),b(24,n3),c(n1,n3)
2724c
2725 do j=1,n3
2726 do i=1,n1
2727 c(i,j) = a(i,1)*b(1,j)
2728 $ + a(i,2)*b(2,j)
2729 $ + a(i,3)*b(3,j)
2730 $ + a(i,4)*b(4,j)
2731 $ + a(i,5)*b(5,j)
2732 $ + a(i,6)*b(6,j)
2733 $ + a(i,7)*b(7,j)
2734 $ + a(i,8)*b(8,j)
2735 $ + a(i,9)*b(9,j)
2736 $ + a(i,10)*b(10,j)
2737 $ + a(i,11)*b(11,j)
2738 $ + a(i,12)*b(12,j)
2739 $ + a(i,13)*b(13,j)
2740 $ + a(i,14)*b(14,j)
2741 $ + a(i,15)*b(15,j)
2742 $ + a(i,16)*b(16,j)
2743 $ + a(i,17)*b(17,j)
2744 $ + a(i,18)*b(18,j)
2745 $ + a(i,19)*b(19,j)
2746 $ + a(i,20)*b(20,j)
2747 $ + a(i,21)*b(21,j)
2748 $ + a(i,22)*b(22,j)
2749 $ + a(i,23)*b(23,j)
2750 $ + a(i,24)*b(24,j)
2751 enddo
2752 enddo
2753 return
2754 end
2755c-----------------------------------------------------------------------
2756 subroutine mxmf3(a,n1,b,n2,c,n3)
2757C-----------------------------------------------------------------------
2758C
2759C Matrix-vector product routine.
2760C NOTE: Use assembly coded routine if available.
2761C
2762C----------------------------------------------------------------------
2763 REAL A(N1,N2),B(N2,N3),C(N1,N3)
2764C
2765 integer wdsize
2766 save wdsize
2767 data wdsize/0/
2768c
2769c First call: determine word size for dgemm/sgemm discrimination, below.
2770c
2771 if (wdsize.eq.0) then
2772 one = 1.0
2773 eps = 1.e-12
2774 wdsize = 8
2775 if (one+eps.eq.1.0) wdsize = 4
2776 endif
2777c
2778 if (n2.le.8) then
2779 if (n2.eq.1) then
2780 call mxmf3_1(a,n1,b,n2,c,n3)
2781 elseif (n2.eq.2) then
2782 call mxmf3_2(a,n1,b,n2,c,n3)
2783 elseif (n2.eq.3) then
2784 call mxmf3_3(a,n1,b,n2,c,n3)
2785 elseif (n2.eq.4) then
2786 call mxmf3_4(a,n1,b,n2,c,n3)
2787 elseif (n2.eq.5) then
2788 call mxmf3_5(a,n1,b,n2,c,n3)
2789 elseif (n2.eq.6) then
2790 call mxmf3_6(a,n1,b,n2,c,n3)
2791 elseif (n2.eq.7) then
2792 call mxmf3_7(a,n1,b,n2,c,n3)
2793 else
2794 call mxmf3_8(a,n1,b,n2,c,n3)
2795 endif
2796 elseif (n2.le.16) then
2797 if (n2.eq.9) then
2798 call mxmf3_9(a,n1,b,n2,c,n3)
2799 elseif (n2.eq.10) then
2800 call mxmf3_10(a,n1,b,n2,c,n3)
2801 elseif (n2.eq.11) then
2802 call mxmf3_11(a,n1,b,n2,c,n3)
2803 elseif (n2.eq.12) then
2804 call mxmf3_12(a,n1,b,n2,c,n3)
2805 elseif (n2.eq.13) then
2806 call mxmf3_13(a,n1,b,n2,c,n3)
2807 elseif (n2.eq.14) then
2808 call mxmf3_14(a,n1,b,n2,c,n3)
2809 elseif (n2.eq.15) then
2810 call mxmf3_15(a,n1,b,n2,c,n3)
2811 else
2812 call mxmf3_16(a,n1,b,n2,c,n3)
2813 endif
2814 elseif (n2.le.24) then
2815 if (n2.eq.17) then
2816 call mxmf3_17(a,n1,b,n2,c,n3)
2817 elseif (n2.eq.18) then
2818 call mxmf3_18(a,n1,b,n2,c,n3)
2819 elseif (n2.eq.19) then
2820 call mxmf3_19(a,n1,b,n2,c,n3)
2821 elseif (n2.eq.20) then
2822 call mxmf3_20(a,n1,b,n2,c,n3)
2823 elseif (n2.eq.21) then
2824 call mxmf3_21(a,n1,b,n2,c,n3)
2825 elseif (n2.eq.22) then
2826 call mxmf3_22(a,n1,b,n2,c,n3)
2827 elseif (n2.eq.23) then
2828 call mxmf3_23(a,n1,b,n2,c,n3)
2829 elseif (n2.eq.24) then
2830 call mxmf3_24(a,n1,b,n2,c,n3)
2831 endif
2832 else
2833c
2834 one=1.0
2835 zero=0.0
2836 call dgemm( 'N','N',n1,n3,n2,ONE,A,N1,B,N2,ZERO,C,N1)
2837c
2838c N0=N1*N3
2839c DO 10 I=1,N0
2840c C(I,1)=0.
2841c 10 CONTINUE
2842c DO 100 J=1,N3
2843c DO 100 K=1,N2
2844c BB=B(K,J)
2845c DO 100 I=1,N1
2846c C(I,J)=C(I,J)+A(I,K)*BB
2847c 100 CONTINUE
2848
2849 endif
2850 return
2851 end
2852c-----------------------------------------------------------------------
2853 subroutine mxmf3_1(a,n1,b,n2,c,n3)
2854c
2855 real a(n1,1),b(1,n3),c(n1,n3)
2856c
2857 do i=1,n1
2858 do j=1,n3
2859 c(i,j) = a(i,1)*b(1,j)
2860 enddo
2861 enddo
2862 return
2863 end
2864c-----------------------------------------------------------------------
2865 subroutine mxmf3_2(a,n1,b,n2,c,n3)
2866c
2867 real a(n1,2),b(2,n3),c(n1,n3)
2868c
2869 do i=1,n1
2870 do j=1,n3
2871 c(i,j) = a(i,1)*b(1,j)
2872 $ + a(i,2)*b(2,j)
2873 enddo
2874 enddo
2875 return
2876 end
2877c-----------------------------------------------------------------------
2878 subroutine mxmf3_3(a,n1,b,n2,c,n3)
2879c
2880 real a(n1,3),b(3,n3),c(n1,n3)
2881c
2882 do i=1,n1
2883 do j=1,n3
2884 c(i,j) = a(i,1)*b(1,j)
2885 $ + a(i,2)*b(2,j)
2886 $ + a(i,3)*b(3,j)
2887 enddo
2888 enddo
2889 return
2890 end
2891c-----------------------------------------------------------------------
2892 subroutine mxmf3_4(a,n1,b,n2,c,n3)
2893c
2894 real a(n1,4),b(4,n3),c(n1,n3)
2895c
2896 do i=1,n1
2897 do j=1,n3
2898 c(i,j) = a(i,1)*b(1,j)
2899 $ + a(i,2)*b(2,j)
2900 $ + a(i,3)*b(3,j)
2901 $ + a(i,4)*b(4,j)
2902 enddo
2903 enddo
2904 return
2905 end
2906c-----------------------------------------------------------------------
2907 subroutine mxmf3_5(a,n1,b,n2,c,n3)
2908c
2909 real a(n1,5),b(5,n3),c(n1,n3)
2910c
2911 do i=1,n1
2912 do j=1,n3
2913 c(i,j) = a(i,1)*b(1,j)
2914 $ + a(i,2)*b(2,j)
2915 $ + a(i,3)*b(3,j)
2916 $ + a(i,4)*b(4,j)
2917 $ + a(i,5)*b(5,j)
2918 enddo
2919 enddo
2920 return
2921 end
2922c-----------------------------------------------------------------------
2923 subroutine mxmf3_6(a,n1,b,n2,c,n3)
2924c
2925 real a(n1,6),b(6,n3),c(n1,n3)
2926c
2927 do i=1,n1
2928 do j=1,n3
2929 c(i,j) = a(i,1)*b(1,j)
2930 $ + a(i,2)*b(2,j)
2931 $ + a(i,3)*b(3,j)
2932 $ + a(i,4)*b(4,j)
2933 $ + a(i,5)*b(5,j)
2934 $ + a(i,6)*b(6,j)
2935 enddo
2936 enddo
2937 return
2938 end
2939c-----------------------------------------------------------------------
2940 subroutine mxmf3_7(a,n1,b,n2,c,n3)
2941c
2942 real a(n1,7),b(7,n3),c(n1,n3)
2943c
2944 do i=1,n1
2945 do j=1,n3
2946 c(i,j) = a(i,1)*b(1,j)
2947 $ + a(i,2)*b(2,j)
2948 $ + a(i,3)*b(3,j)
2949 $ + a(i,4)*b(4,j)
2950 $ + a(i,5)*b(5,j)
2951 $ + a(i,6)*b(6,j)
2952 $ + a(i,7)*b(7,j)
2953 enddo
2954 enddo
2955 return
2956 end
2957c-----------------------------------------------------------------------
2958 subroutine mxmf3_8(a,n1,b,n2,c,n3)
2959c
2960 real a(n1,8),b(8,n3),c(n1,n3)
2961c
2962 do i=1,n1
2963 do j=1,n3
2964 c(i,j) = a(i,1)*b(1,j)
2965 $ + a(i,2)*b(2,j)
2966 $ + a(i,3)*b(3,j)
2967 $ + a(i,4)*b(4,j)
2968 $ + a(i,5)*b(5,j)
2969 $ + a(i,6)*b(6,j)
2970 $ + a(i,7)*b(7,j)
2971 $ + a(i,8)*b(8,j)
2972 enddo
2973 enddo
2974 return
2975 end
2976c-----------------------------------------------------------------------
2977 subroutine mxmf3_9(a,n1,b,n2,c,n3)
2978c
2979 real a(n1,9),b(9,n3),c(n1,n3)
2980c
2981 do i=1,n1
2982 do j=1,n3
2983 c(i,j) = a(i,1)*b(1,j)
2984 $ + a(i,2)*b(2,j)
2985 $ + a(i,3)*b(3,j)
2986 $ + a(i,4)*b(4,j)
2987 $ + a(i,5)*b(5,j)
2988 $ + a(i,6)*b(6,j)
2989 $ + a(i,7)*b(7,j)
2990 $ + a(i,8)*b(8,j)
2991 $ + a(i,9)*b(9,j)
2992 enddo
2993 enddo
2994 return
2995 end
2996c-----------------------------------------------------------------------
2997 subroutine mxmf3_10(a,n1,b,n2,c,n3)
2998c
2999 real a(n1,10),b(10,n3),c(n1,n3)
3000c
3001 do i=1,n1
3002 do j=1,n3
3003 c(i,j) = a(i,1)*b(1,j)
3004 $ + a(i,2)*b(2,j)
3005 $ + a(i,3)*b(3,j)
3006 $ + a(i,4)*b(4,j)
3007 $ + a(i,5)*b(5,j)
3008 $ + a(i,6)*b(6,j)
3009 $ + a(i,7)*b(7,j)
3010 $ + a(i,8)*b(8,j)
3011 $ + a(i,9)*b(9,j)
3012 $ + a(i,10)*b(10,j)
3013 enddo
3014 enddo
3015 return
3016 end
3017c-----------------------------------------------------------------------
3018 subroutine mxmf3_11(a,n1,b,n2,c,n3)
3019c
3020 real a(n1,11),b(11,n3),c(n1,n3)
3021c
3022 do i=1,n1
3023 do j=1,n3
3024 c(i,j) = a(i,1)*b(1,j)
3025 $ + a(i,2)*b(2,j)
3026 $ + a(i,3)*b(3,j)
3027 $ + a(i,4)*b(4,j)
3028 $ + a(i,5)*b(5,j)
3029 $ + a(i,6)*b(6,j)
3030 $ + a(i,7)*b(7,j)
3031 $ + a(i,8)*b(8,j)
3032 $ + a(i,9)*b(9,j)
3033 $ + a(i,10)*b(10,j)
3034 $ + a(i,11)*b(11,j)
3035 enddo
3036 enddo
3037 return
3038 end
3039c-----------------------------------------------------------------------
3040 subroutine mxmf3_12(a,n1,b,n2,c,n3)
3041c
3042 real a(n1,12),b(12,n3),c(n1,n3)
3043c
3044 do i=1,n1
3045 do j=1,n3
3046 c(i,j) = a(i,1)*b(1,j)
3047 $ + a(i,2)*b(2,j)
3048 $ + a(i,3)*b(3,j)
3049 $ + a(i,4)*b(4,j)
3050 $ + a(i,5)*b(5,j)
3051 $ + a(i,6)*b(6,j)
3052 $ + a(i,7)*b(7,j)
3053 $ + a(i,8)*b(8,j)
3054 $ + a(i,9)*b(9,j)
3055 $ + a(i,10)*b(10,j)
3056 $ + a(i,11)*b(11,j)
3057 $ + a(i,12)*b(12,j)
3058 enddo
3059 enddo
3060 return
3061 end
3062c-----------------------------------------------------------------------
3063 subroutine mxmf3_13(a,n1,b,n2,c,n3)
3064c
3065 real a(n1,13),b(13,n3),c(n1,n3)
3066c
3067 do i=1,n1
3068 do j=1,n3
3069 c(i,j) = a(i,1)*b(1,j)
3070 $ + a(i,2)*b(2,j)
3071 $ + a(i,3)*b(3,j)
3072 $ + a(i,4)*b(4,j)
3073 $ + a(i,5)*b(5,j)
3074 $ + a(i,6)*b(6,j)
3075 $ + a(i,7)*b(7,j)
3076 $ + a(i,8)*b(8,j)
3077 $ + a(i,9)*b(9,j)
3078 $ + a(i,10)*b(10,j)
3079 $ + a(i,11)*b(11,j)
3080 $ + a(i,12)*b(12,j)
3081 $ + a(i,13)*b(13,j)
3082 enddo
3083 enddo
3084 return
3085 end
3086c-----------------------------------------------------------------------
3087 subroutine mxmf3_14(a,n1,b,n2,c,n3)
3088c
3089 real a(n1,14),b(14,n3),c(n1,n3)
3090c
3091 do i=1,n1
3092 do j=1,n3
3093 c(i,j) = a(i,1)*b(1,j)
3094 $ + a(i,2)*b(2,j)
3095 $ + a(i,3)*b(3,j)
3096 $ + a(i,4)*b(4,j)
3097 $ + a(i,5)*b(5,j)
3098 $ + a(i,6)*b(6,j)
3099 $ + a(i,7)*b(7,j)
3100 $ + a(i,8)*b(8,j)
3101 $ + a(i,9)*b(9,j)
3102 $ + a(i,10)*b(10,j)
3103 $ + a(i,11)*b(11,j)
3104 $ + a(i,12)*b(12,j)
3105 $ + a(i,13)*b(13,j)
3106 $ + a(i,14)*b(14,j)
3107 enddo
3108 enddo
3109 return
3110 end
3111c-----------------------------------------------------------------------
3112 subroutine mxmf3_15(a,n1,b,n2,c,n3)
3113c
3114 real a(n1,15),b(15,n3),c(n1,n3)
3115c
3116 do i=1,n1
3117 do j=1,n3
3118 c(i,j) = a(i,1)*b(1,j)
3119 $ + a(i,2)*b(2,j)
3120 $ + a(i,3)*b(3,j)
3121 $ + a(i,4)*b(4,j)
3122 $ + a(i,5)*b(5,j)
3123 $ + a(i,6)*b(6,j)
3124 $ + a(i,7)*b(7,j)
3125 $ + a(i,8)*b(8,j)
3126 $ + a(i,9)*b(9,j)
3127 $ + a(i,10)*b(10,j)
3128 $ + a(i,11)*b(11,j)
3129 $ + a(i,12)*b(12,j)
3130 $ + a(i,13)*b(13,j)
3131 $ + a(i,14)*b(14,j)
3132 $ + a(i,15)*b(15,j)
3133 enddo
3134 enddo
3135 return
3136 end
3137c-----------------------------------------------------------------------
3138 subroutine mxmf3_16(a,n1,b,n2,c,n3)
3139c
3140 real a(n1,16),b(16,n3),c(n1,n3)
3141c
3142 do i=1,n1
3143 do j=1,n3
3144 c(i,j) = a(i,1)*b(1,j)
3145 $ + a(i,2)*b(2,j)
3146 $ + a(i,3)*b(3,j)
3147 $ + a(i,4)*b(4,j)
3148 $ + a(i,5)*b(5,j)
3149 $ + a(i,6)*b(6,j)
3150 $ + a(i,7)*b(7,j)
3151 $ + a(i,8)*b(8,j)
3152 $ + a(i,9)*b(9,j)
3153 $ + a(i,10)*b(10,j)
3154 $ + a(i,11)*b(11,j)
3155 $ + a(i,12)*b(12,j)
3156 $ + a(i,13)*b(13,j)
3157 $ + a(i,14)*b(14,j)
3158 $ + a(i,15)*b(15,j)
3159 $ + a(i,16)*b(16,j)
3160 enddo
3161 enddo
3162 return
3163 end
3164c-----------------------------------------------------------------------
3165 subroutine mxmf3_17(a,n1,b,n2,c,n3)
3166c
3167 real a(n1,17),b(17,n3),c(n1,n3)
3168c
3169 do i=1,n1
3170 do j=1,n3
3171 c(i,j) = a(i,1)*b(1,j)
3172 $ + a(i,2)*b(2,j)
3173 $ + a(i,3)*b(3,j)
3174 $ + a(i,4)*b(4,j)
3175 $ + a(i,5)*b(5,j)
3176 $ + a(i,6)*b(6,j)
3177 $ + a(i,7)*b(7,j)
3178 $ + a(i,8)*b(8,j)
3179 $ + a(i,9)*b(9,j)
3180 $ + a(i,10)*b(10,j)
3181 $ + a(i,11)*b(11,j)
3182 $ + a(i,12)*b(12,j)
3183 $ + a(i,13)*b(13,j)
3184 $ + a(i,14)*b(14,j)
3185 $ + a(i,15)*b(15,j)
3186 $ + a(i,16)*b(16,j)
3187 $ + a(i,17)*b(17,j)
3188 enddo
3189 enddo
3190 return
3191 end
3192c-----------------------------------------------------------------------
3193 subroutine mxmf3_18(a,n1,b,n2,c,n3)
3194c
3195 real a(n1,18),b(18,n3),c(n1,n3)
3196c
3197 do i=1,n1
3198 do j=1,n3
3199 c(i,j) = a(i,1)*b(1,j)
3200 $ + a(i,2)*b(2,j)
3201 $ + a(i,3)*b(3,j)
3202 $ + a(i,4)*b(4,j)
3203 $ + a(i,5)*b(5,j)
3204 $ + a(i,6)*b(6,j)
3205 $ + a(i,7)*b(7,j)
3206 $ + a(i,8)*b(8,j)
3207 $ + a(i,9)*b(9,j)
3208 $ + a(i,10)*b(10,j)
3209 $ + a(i,11)*b(11,j)
3210 $ + a(i,12)*b(12,j)
3211 $ + a(i,13)*b(13,j)
3212 $ + a(i,14)*b(14,j)
3213 $ + a(i,15)*b(15,j)
3214 $ + a(i,16)*b(16,j)
3215 $ + a(i,17)*b(17,j)
3216 $ + a(i,18)*b(18,j)
3217 enddo
3218 enddo
3219 return
3220 end
3221c-----------------------------------------------------------------------
3222 subroutine mxmf3_19(a,n1,b,n2,c,n3)
3223c
3224 real a(n1,19),b(19,n3),c(n1,n3)
3225c
3226 do i=1,n1
3227 do j=1,n3
3228 c(i,j) = a(i,1)*b(1,j)
3229 $ + a(i,2)*b(2,j)
3230 $ + a(i,3)*b(3,j)
3231 $ + a(i,4)*b(4,j)
3232 $ + a(i,5)*b(5,j)
3233 $ + a(i,6)*b(6,j)
3234 $ + a(i,7)*b(7,j)
3235 $ + a(i,8)*b(8,j)
3236 $ + a(i,9)*b(9,j)
3237 $ + a(i,10)*b(10,j)
3238 $ + a(i,11)*b(11,j)
3239 $ + a(i,12)*b(12,j)
3240 $ + a(i,13)*b(13,j)
3241 $ + a(i,14)*b(14,j)
3242 $ + a(i,15)*b(15,j)
3243 $ + a(i,16)*b(16,j)
3244 $ + a(i,17)*b(17,j)
3245 $ + a(i,18)*b(18,j)
3246 $ + a(i,19)*b(19,j)
3247 enddo
3248 enddo
3249 return
3250 end
3251c-----------------------------------------------------------------------
3252 subroutine mxmf3_20(a,n1,b,n2,c,n3)
3253c
3254 real a(n1,20),b(20,n3),c(n1,n3)
3255c
3256 do i=1,n1
3257 do j=1,n3
3258 c(i,j) = a(i,1)*b(1,j)
3259 $ + a(i,2)*b(2,j)
3260 $ + a(i,3)*b(3,j)
3261 $ + a(i,4)*b(4,j)
3262 $ + a(i,5)*b(5,j)
3263 $ + a(i,6)*b(6,j)
3264 $ + a(i,7)*b(7,j)
3265 $ + a(i,8)*b(8,j)
3266 $ + a(i,9)*b(9,j)
3267 $ + a(i,10)*b(10,j)
3268 $ + a(i,11)*b(11,j)
3269 $ + a(i,12)*b(12,j)
3270 $ + a(i,13)*b(13,j)
3271 $ + a(i,14)*b(14,j)
3272 $ + a(i,15)*b(15,j)
3273 $ + a(i,16)*b(16,j)
3274 $ + a(i,17)*b(17,j)
3275 $ + a(i,18)*b(18,j)
3276 $ + a(i,19)*b(19,j)
3277 $ + a(i,20)*b(20,j)
3278 enddo
3279 enddo
3280 return
3281 end
3282c-----------------------------------------------------------------------
3283 subroutine mxmf3_21(a,n1,b,n2,c,n3)
3284c
3285 real a(n1,21),b(21,n3),c(n1,n3)
3286c
3287 do i=1,n1
3288 do j=1,n3
3289 c(i,j) = a(i,1)*b(1,j)
3290 $ + a(i,2)*b(2,j)
3291 $ + a(i,3)*b(3,j)
3292 $ + a(i,4)*b(4,j)
3293 $ + a(i,5)*b(5,j)
3294 $ + a(i,6)*b(6,j)
3295 $ + a(i,7)*b(7,j)
3296 $ + a(i,8)*b(8,j)
3297 $ + a(i,9)*b(9,j)
3298 $ + a(i,10)*b(10,j)
3299 $ + a(i,11)*b(11,j)
3300 $ + a(i,12)*b(12,j)
3301 $ + a(i,13)*b(13,j)
3302 $ + a(i,14)*b(14,j)
3303 $ + a(i,15)*b(15,j)
3304 $ + a(i,16)*b(16,j)
3305 $ + a(i,17)*b(17,j)
3306 $ + a(i,18)*b(18,j)
3307 $ + a(i,19)*b(19,j)
3308 $ + a(i,20)*b(20,j)
3309 $ + a(i,21)*b(21,j)
3310 enddo
3311 enddo
3312 return
3313 end
3314c-----------------------------------------------------------------------
3315 subroutine mxmf3_22(a,n1,b,n2,c,n3)
3316c
3317 real a(n1,22),b(22,n3),c(n1,n3)
3318c
3319 do i=1,n1
3320 do j=1,n3
3321 c(i,j) = a(i,1)*b(1,j)
3322 $ + a(i,2)*b(2,j)
3323 $ + a(i,3)*b(3,j)
3324 $ + a(i,4)*b(4,j)
3325 $ + a(i,5)*b(5,j)
3326 $ + a(i,6)*b(6,j)
3327 $ + a(i,7)*b(7,j)
3328 $ + a(i,8)*b(8,j)
3329 $ + a(i,9)*b(9,j)
3330 $ + a(i,10)*b(10,j)
3331 $ + a(i,11)*b(11,j)
3332 $ + a(i,12)*b(12,j)
3333 $ + a(i,13)*b(13,j)
3334 $ + a(i,14)*b(14,j)
3335 $ + a(i,15)*b(15,j)
3336 $ + a(i,16)*b(16,j)
3337 $ + a(i,17)*b(17,j)
3338 $ + a(i,18)*b(18,j)
3339 $ + a(i,19)*b(19,j)
3340 $ + a(i,20)*b(20,j)
3341 $ + a(i,21)*b(21,j)
3342 $ + a(i,22)*b(22,j)
3343 enddo
3344 enddo
3345 return
3346 end
3347c-----------------------------------------------------------------------
3348 subroutine mxmf3_23(a,n1,b,n2,c,n3)
3349c
3350 real a(n1,23),b(23,n3),c(n1,n3)
3351c
3352 do i=1,n1
3353 do j=1,n3
3354 c(i,j) = a(i,1)*b(1,j)
3355 $ + a(i,2)*b(2,j)
3356 $ + a(i,3)*b(3,j)
3357 $ + a(i,4)*b(4,j)
3358 $ + a(i,5)*b(5,j)
3359 $ + a(i,6)*b(6,j)
3360 $ + a(i,7)*b(7,j)
3361 $ + a(i,8)*b(8,j)
3362 $ + a(i,9)*b(9,j)
3363 $ + a(i,10)*b(10,j)
3364 $ + a(i,11)*b(11,j)
3365 $ + a(i,12)*b(12,j)
3366 $ + a(i,13)*b(13,j)
3367 $ + a(i,14)*b(14,j)
3368 $ + a(i,15)*b(15,j)
3369 $ + a(i,16)*b(16,j)
3370 $ + a(i,17)*b(17,j)
3371 $ + a(i,18)*b(18,j)
3372 $ + a(i,19)*b(19,j)
3373 $ + a(i,20)*b(20,j)
3374 $ + a(i,21)*b(21,j)
3375 $ + a(i,22)*b(22,j)
3376 $ + a(i,23)*b(23,j)
3377 enddo
3378 enddo
3379 return
3380 end
3381c-----------------------------------------------------------------------
3382 subroutine mxmf3_24(a,n1,b,n2,c,n3)
3383c
3384 real a(n1,24),b(24,n3),c(n1,n3)
3385c
3386 do i=1,n1
3387 do j=1,n3
3388 c(i,j) = a(i,1)*b(1,j)
3389 $ + a(i,2)*b(2,j)
3390 $ + a(i,3)*b(3,j)
3391 $ + a(i,4)*b(4,j)
3392 $ + a(i,5)*b(5,j)
3393 $ + a(i,6)*b(6,j)
3394 $ + a(i,7)*b(7,j)
3395 $ + a(i,8)*b(8,j)
3396 $ + a(i,9)*b(9,j)
3397 $ + a(i,10)*b(10,j)
3398 $ + a(i,11)*b(11,j)
3399 $ + a(i,12)*b(12,j)
3400 $ + a(i,13)*b(13,j)
3401 $ + a(i,14)*b(14,j)
3402 $ + a(i,15)*b(15,j)
3403 $ + a(i,16)*b(16,j)
3404 $ + a(i,17)*b(17,j)
3405 $ + a(i,18)*b(18,j)
3406 $ + a(i,19)*b(19,j)
3407 $ + a(i,20)*b(20,j)
3408 $ + a(i,21)*b(21,j)
3409 $ + a(i,22)*b(22,j)
3410 $ + a(i,23)*b(23,j)
3411 $ + a(i,24)*b(24,j)
3412 enddo
3413 enddo
3414 return
3415 end
3416c-----------------------------------------------------------------------
3417 subroutine mxm44(a,n1,b,n2,c,n3)
3418C-----------------------------------------------------------------------
3419C
3420C NOTE -- this code has been set up with the "mxmf3" routine
3421c referenced in memtime.f. On most machines, the f2
3422c and f3 versions give the same performance (f2 is the
3423c nekton standard). On the t3e, f3 is noticeably faster.
3424c pff 10/5/98
3425C
3426C
3427C Matrix-vector product routine.
3428C NOTE: Use assembly coded routine if available.
3429C
3430C----------------------------------------------------------------------
3431 REAL A(N1,N2),B(N2,N3),C(N1,N3)
3432c
3433 if (n2.eq.1) then
3434 call mxm44_2_t(a,n1,b,2,c,n3)
3435 elseif (n2.eq.2) then
3436 call mxm44_2_t(a,n1,b,n2,c,n3)
3437 else
3438 call mxm44_0_t(a,n1,b,n2,c,n3)
3439 endif
3440c
3441 return
3442 end
3443c
3444c-----------------------------------------------------------------------
3445 subroutine mxm44_0_t(a, m, b, k, c, n)
3446* subroutine matmul44(m, n, k, a, lda, b, ldb, c, ldc)
3447* real*8 a(lda,k), b(ldb,n), c(ldc,n)
3448 real a(m,k), b(k,n), c(m,n)
3449 real s11, s12, s13, s14, s21, s22, s23, s24
3450 real s31, s32, s33, s34, s41, s42, s43, s44
3451c
3452c matrix multiply with a 4x4 pencil
3453c
3454
3455 mresid = iand(m,3)
3456 nresid = iand(n,3)
3457 m1 = m - mresid + 1
3458 n1 = n - nresid + 1
3459
3460 do i=1,m-mresid,4
3461 do j=1,n-nresid,4
3462 s11 = 0.0d0
3463 s21 = 0.0d0
3464 s31 = 0.0d0
3465 s41 = 0.0d0
3466 s12 = 0.0d0
3467 s22 = 0.0d0
3468 s32 = 0.0d0
3469 s42 = 0.0d0
3470 s13 = 0.0d0
3471 s23 = 0.0d0
3472 s33 = 0.0d0
3473 s43 = 0.0d0
3474 s14 = 0.0d0
3475 s24 = 0.0d0
3476 s34 = 0.0d0
3477 s44 = 0.0d0
3478 do l=1,k
3479 s11 = s11 + a(i,l)*b(l,j)
3480 s12 = s12 + a(i,l)*b(l,j+1)
3481 s13 = s13 + a(i,l)*b(l,j+2)
3482 s14 = s14 + a(i,l)*b(l,j+3)
3483
3484 s21 = s21 + a(i+1,l)*b(l,j)
3485 s22 = s22 + a(i+1,l)*b(l,j+1)
3486 s23 = s23 + a(i+1,l)*b(l,j+2)
3487 s24 = s24 + a(i+1,l)*b(l,j+3)
3488
3489 s31 = s31 + a(i+2,l)*b(l,j)
3490 s32 = s32 + a(i+2,l)*b(l,j+1)
3491 s33 = s33 + a(i+2,l)*b(l,j+2)
3492 s34 = s34 + a(i+2,l)*b(l,j+3)
3493
3494 s41 = s41 + a(i+3,l)*b(l,j)
3495 s42 = s42 + a(i+3,l)*b(l,j+1)
3496 s43 = s43 + a(i+3,l)*b(l,j+2)
3497 s44 = s44 + a(i+3,l)*b(l,j+3)
3498 enddo
3499 c(i,j) = s11
3500 c(i,j+1) = s12
3501 c(i,j+2) = s13
3502 c(i,j+3) = s14
3503
3504 c(i+1,j) = s21
3505 c(i+2,j) = s31
3506 c(i+3,j) = s41
3507
3508 c(i+1,j+1) = s22
3509 c(i+2,j+1) = s32
3510 c(i+3,j+1) = s42
3511
3512 c(i+1,j+2) = s23
3513 c(i+2,j+2) = s33
3514 c(i+3,j+2) = s43
3515
3516 c(i+1,j+3) = s24
3517 c(i+2,j+3) = s34
3518 c(i+3,j+3) = s44
3519 enddo
3520* Residual when n is not multiple of 4
3521 if (nresid .ne. 0) then
3522 if (nresid .eq. 1) then
3523 s11 = 0.0d0
3524 s21 = 0.0d0
3525 s31 = 0.0d0
3526 s41 = 0.0d0
3527 do l=1,k
3528 s11 = s11 + a(i,l)*b(l,n)
3529 s21 = s21 + a(i+1,l)*b(l,n)
3530 s31 = s31 + a(i+2,l)*b(l,n)
3531 s41 = s41 + a(i+3,l)*b(l,n)
3532 enddo
3533 c(i,n) = s11
3534 c(i+1,n) = s21
3535 c(i+2,n) = s31
3536 c(i+3,n) = s41
3537 elseif (nresid .eq. 2) then
3538 s11 = 0.0d0
3539 s21 = 0.0d0
3540 s31 = 0.0d0
3541 s41 = 0.0d0
3542 s12 = 0.0d0
3543 s22 = 0.0d0
3544 s32 = 0.0d0
3545 s42 = 0.0d0
3546 do l=1,k
3547 s11 = s11 + a(i,l)*b(l,j)
3548 s12 = s12 + a(i,l)*b(l,j+1)
3549
3550 s21 = s21 + a(i+1,l)*b(l,j)
3551 s22 = s22 + a(i+1,l)*b(l,j+1)
3552
3553 s31 = s31 + a(i+2,l)*b(l,j)
3554 s32 = s32 + a(i+2,l)*b(l,j+1)
3555
3556 s41 = s41 + a(i+3,l)*b(l,j)
3557 s42 = s42 + a(i+3,l)*b(l,j+1)
3558 enddo
3559 c(i,j) = s11
3560 c(i,j+1) = s12
3561
3562 c(i+1,j) = s21
3563 c(i+2,j) = s31
3564 c(i+3,j) = s41
3565
3566 c(i+1,j+1) = s22
3567 c(i+2,j+1) = s32
3568 c(i+3,j+1) = s42
3569 else
3570 s11 = 0.0d0
3571 s21 = 0.0d0
3572 s31 = 0.0d0
3573 s41 = 0.0d0
3574 s12 = 0.0d0
3575 s22 = 0.0d0
3576 s32 = 0.0d0
3577 s42 = 0.0d0
3578 s13 = 0.0d0
3579 s23 = 0.0d0
3580 s33 = 0.0d0
3581 s43 = 0.0d0
3582 do l=1,k
3583 s11 = s11 + a(i,l)*b(l,j)
3584 s12 = s12 + a(i,l)*b(l,j+1)
3585 s13 = s13 + a(i,l)*b(l,j+2)
3586
3587 s21 = s21 + a(i+1,l)*b(l,j)
3588 s22 = s22 + a(i+1,l)*b(l,j+1)
3589 s23 = s23 + a(i+1,l)*b(l,j+2)
3590
3591 s31 = s31 + a(i+2,l)*b(l,j)
3592 s32 = s32 + a(i+2,l)*b(l,j+1)
3593 s33 = s33 + a(i+2,l)*b(l,j+2)
3594
3595 s41 = s41 + a(i+3,l)*b(l,j)
3596 s42 = s42 + a(i+3,l)*b(l,j+1)
3597 s43 = s43 + a(i+3,l)*b(l,j+2)
3598 enddo
3599 c(i,j) = s11
3600 c(i+1,j) = s21
3601 c(i+2,j) = s31
3602 c(i+3,j) = s41
3603 c(i,j+1) = s12
3604 c(i+1,j+1) = s22
3605 c(i+2,j+1) = s32
3606 c(i+3,j+1) = s42
3607 c(i,j+2) = s13
3608 c(i+1,j+2) = s23
3609 c(i+2,j+2) = s33
3610 c(i+3,j+2) = s43
3611 endif
3612 endif
3613 enddo
3614
3615* Residual when m is not multiple of 4
3616 if (mresid .eq. 0) then
3617 return
3618 elseif (mresid .eq. 1) then
3619 do j=1,n-nresid,4
3620 s11 = 0.0d0
3621 s12 = 0.0d0
3622 s13 = 0.0d0
3623 s14 = 0.0d0
3624 do l=1,k
3625 s11 = s11 + a(m,l)*b(l,j)
3626 s12 = s12 + a(m,l)*b(l,j+1)
3627 s13 = s13 + a(m,l)*b(l,j+2)
3628 s14 = s14 + a(m,l)*b(l,j+3)
3629 enddo
3630 c(m,j) = s11
3631 c(m,j+1) = s12
3632 c(m,j+2) = s13
3633 c(m,j+3) = s14
3634 enddo
3635* mresid is 1, check nresid
3636 if (nresid .eq. 0) then
3637 return
3638 elseif (nresid .eq. 1) then
3639 s11 = 0.0d0
3640 do l=1,k
3641 s11 = s11 + a(m,l)*b(l,n)
3642 enddo
3643 c(m,n) = s11
3644 return
3645 elseif (nresid .eq. 2) then
3646 s11 = 0.0d0
3647 s12 = 0.0d0
3648 do l=1,k
3649 s11 = s11 + a(m,l)*b(l,n-1)
3650 s12 = s12 + a(m,l)*b(l,n)
3651 enddo
3652 c(m,n-1) = s11
3653 c(m,n) = s12
3654 return
3655 else
3656 s11 = 0.0d0
3657 s12 = 0.0d0
3658 s13 = 0.0d0
3659 do l=1,k
3660 s11 = s11 + a(m,l)*b(l,n-2)
3661 s12 = s12 + a(m,l)*b(l,n-1)
3662 s13 = s13 + a(m,l)*b(l,n)
3663 enddo
3664 c(m,n-2) = s11
3665 c(m,n-1) = s12
3666 c(m,n) = s13
3667 return
3668 endif
3669 elseif (mresid .eq. 2) then
3670 do j=1,n-nresid,4
3671 s11 = 0.0d0
3672 s12 = 0.0d0
3673 s13 = 0.0d0
3674 s14 = 0.0d0
3675 s21 = 0.0d0
3676 s22 = 0.0d0
3677 s23 = 0.0d0
3678 s24 = 0.0d0
3679 do l=1,k
3680 s11 = s11 + a(m-1,l)*b(l,j)
3681 s12 = s12 + a(m-1,l)*b(l,j+1)
3682 s13 = s13 + a(m-1,l)*b(l,j+2)
3683 s14 = s14 + a(m-1,l)*b(l,j+3)
3684
3685 s21 = s21 + a(m,l)*b(l,j)
3686 s22 = s22 + a(m,l)*b(l,j+1)
3687 s23 = s23 + a(m,l)*b(l,j+2)
3688 s24 = s24 + a(m,l)*b(l,j+3)
3689 enddo
3690 c(m-1,j) = s11
3691 c(m-1,j+1) = s12
3692 c(m-1,j+2) = s13
3693 c(m-1,j+3) = s14
3694 c(m,j) = s21
3695 c(m,j+1) = s22
3696 c(m,j+2) = s23
3697 c(m,j+3) = s24
3698 enddo
3699* mresid is 2, check nresid
3700 if (nresid .eq. 0) then
3701 return
3702 elseif (nresid .eq. 1) then
3703 s11 = 0.0d0
3704 s21 = 0.0d0
3705 do l=1,k
3706 s11 = s11 + a(m-1,l)*b(l,n)
3707 s21 = s21 + a(m,l)*b(l,n)
3708 enddo
3709 c(m-1,n) = s11
3710 c(m,n) = s21
3711 return
3712 elseif (nresid .eq. 2) then
3713 s11 = 0.0d0
3714 s21 = 0.0d0
3715 s12 = 0.0d0
3716 s22 = 0.0d0
3717 do l=1,k
3718 s11 = s11 + a(m-1,l)*b(l,n-1)
3719 s12 = s12 + a(m-1,l)*b(l,n)
3720 s21 = s21 + a(m,l)*b(l,n-1)
3721 s22 = s22 + a(m,l)*b(l,n)
3722 enddo
3723 c(m-1,n-1) = s11
3724 c(m-1,n) = s12
3725 c(m,n-1) = s21
3726 c(m,n) = s22
3727 return
3728 else
3729 s11 = 0.0d0
3730 s21 = 0.0d0
3731 s12 = 0.0d0
3732 s22 = 0.0d0
3733 s13 = 0.0d0
3734 s23 = 0.0d0
3735 do l=1,k
3736 s11 = s11 + a(m-1,l)*b(l,n-2)
3737 s12 = s12 + a(m-1,l)*b(l,n-1)
3738 s13 = s13 + a(m-1,l)*b(l,n)
3739 s21 = s21 + a(m,l)*b(l,n-2)
3740 s22 = s22 + a(m,l)*b(l,n-1)
3741 s23 = s23 + a(m,l)*b(l,n)
3742 enddo
3743 c(m-1,n-2) = s11
3744 c(m-1,n-1) = s12
3745 c(m-1,n) = s13
3746 c(m,n-2) = s21
3747 c(m,n-1) = s22
3748 c(m,n) = s23
3749 return
3750 endif
3751 else
3752* mresid is 3
3753 do j=1,n-nresid,4
3754 s11 = 0.0d0
3755 s21 = 0.0d0
3756 s31 = 0.0d0
3757
3758 s12 = 0.0d0
3759 s22 = 0.0d0
3760 s32 = 0.0d0
3761
3762 s13 = 0.0d0
3763 s23 = 0.0d0
3764 s33 = 0.0d0
3765
3766 s14 = 0.0d0
3767 s24 = 0.0d0
3768 s34 = 0.0d0
3769
3770 do l=1,k
3771 s11 = s11 + a(m-2,l)*b(l,j)
3772 s12 = s12 + a(m-2,l)*b(l,j+1)
3773 s13 = s13 + a(m-2,l)*b(l,j+2)
3774 s14 = s14 + a(m-2,l)*b(l,j+3)
3775
3776 s21 = s21 + a(m-1,l)*b(l,j)
3777 s22 = s22 + a(m-1,l)*b(l,j+1)
3778 s23 = s23 + a(m-1,l)*b(l,j+2)
3779 s24 = s24 + a(m-1,l)*b(l,j+3)
3780
3781 s31 = s31 + a(m,l)*b(l,j)
3782 s32 = s32 + a(m,l)*b(l,j+1)
3783 s33 = s33 + a(m,l)*b(l,j+2)
3784 s34 = s34 + a(m,l)*b(l,j+3)
3785 enddo
3786 c(m-2,j) = s11
3787 c(m-2,j+1) = s12
3788 c(m-2,j+2) = s13
3789 c(m-2,j+3) = s14
3790
3791 c(m-1,j) = s21
3792 c(m-1,j+1) = s22
3793 c(m-1,j+2) = s23
3794 c(m-1,j+3) = s24
3795
3796 c(m,j) = s31
3797 c(m,j+1) = s32
3798 c(m,j+2) = s33
3799 c(m,j+3) = s34
3800 enddo
3801* mresid is 3, check nresid
3802 if (nresid .eq. 0) then
3803 return
3804 elseif (nresid .eq. 1) then
3805 s11 = 0.0d0
3806 s21 = 0.0d0
3807 s31 = 0.0d0
3808 do l=1,k
3809 s11 = s11 + a(m-2,l)*b(l,n)
3810 s21 = s21 + a(m-1,l)*b(l,n)
3811 s31 = s31 + a(m,l)*b(l,n)
3812 enddo
3813 c(m-2,n) = s11
3814 c(m-1,n) = s21
3815 c(m,n) = s31
3816 return
3817 elseif (nresid .eq. 2) then
3818 s11 = 0.0d0
3819 s21 = 0.0d0
3820 s31 = 0.0d0
3821 s12 = 0.0d0
3822 s22 = 0.0d0
3823 s32 = 0.0d0
3824 do l=1,k
3825 s11 = s11 + a(m-2,l)*b(l,n-1)
3826 s12 = s12 + a(m-2,l)*b(l,n)
3827 s21 = s21 + a(m-1,l)*b(l,n-1)
3828 s22 = s22 + a(m-1,l)*b(l,n)
3829 s31 = s31 + a(m,l)*b(l,n-1)
3830 s32 = s32 + a(m,l)*b(l,n)
3831 enddo
3832 c(m-2,n-1) = s11
3833 c(m-2,n) = s12
3834 c(m-1,n-1) = s21
3835 c(m-1,n) = s22
3836 c(m,n-1) = s31
3837 c(m,n) = s32
3838 return
3839 else
3840 s11 = 0.0d0
3841 s21 = 0.0d0
3842 s31 = 0.0d0
3843 s12 = 0.0d0
3844 s22 = 0.0d0
3845 s32 = 0.0d0
3846 s13 = 0.0d0
3847 s23 = 0.0d0
3848 s33 = 0.0d0
3849 do l=1,k
3850 s11 = s11 + a(m-2,l)*b(l,n-2)
3851 s12 = s12 + a(m-2,l)*b(l,n-1)
3852 s13 = s13 + a(m-2,l)*b(l,n)
3853 s21 = s21 + a(m-1,l)*b(l,n-2)
3854 s22 = s22 + a(m-1,l)*b(l,n-1)
3855 s23 = s23 + a(m-1,l)*b(l,n)
3856 s31 = s31 + a(m,l)*b(l,n-2)
3857 s32 = s32 + a(m,l)*b(l,n-1)
3858 s33 = s33 + a(m,l)*b(l,n)
3859 enddo
3860 c(m-2,n-2) = s11
3861 c(m-2,n-1) = s12
3862 c(m-2,n) = s13
3863 c(m-1,n-2) = s21
3864 c(m-1,n-1) = s22
3865 c(m-1,n) = s23
3866 c(m,n-2) = s31
3867 c(m,n-1) = s32
3868 c(m,n) = s33
3869 return
3870 endif
3871 endif
3872
3873 return
3874 end
3875c-----------------------------------------------------------------------
3876 subroutine mxm44_2_t(a, m, b, k, c, n)
3877 real a(m,2), b(2,n), c(m,n)
3878
3879 nresid = iand(n,3)
3880 n1 = n - nresid + 1
3881
3882 do j=1,n-nresid,4
3883 do i=1,m
3884 c(i,j) = a(i,1)*b(1,j)
3885 $ + a(i,2)*b(2,j)
3886 c(i,j+1) = a(i,1)*b(1,j+1)
3887 $ + a(i,2)*b(2,j+1)
3888 c(i,j+2) = a(i,1)*b(1,j+2)
3889 $ + a(i,2)*b(2,j+2)
3890 c(i,j+3) = a(i,1)*b(1,j+3)
3891 $ + a(i,2)*b(2,j+3)
3892 enddo
3893 enddo
3894 if (nresid .eq. 0) then
3895 return
3896 elseif (nresid .eq. 1) then
3897 do i=1,m
3898 c(i,n) = a(i,1)*b(1,n)
3899 $ + a(i,2)*b(2,n)
3900 enddo
3901 elseif (nresid .eq. 2) then
3902 do i=1,m
3903 c(i,n-1) = a(i,1)*b(1,n-1)
3904 $ + a(i,2)*b(2,n-1)
3905 c(i,n) = a(i,1)*b(1,n)
3906 $ + a(i,2)*b(2,n)
3907 enddo
3908 else
3909 do i=1,m
3910 c(i,n-2) = a(i,1)*b(1,n-2)
3911 $ + a(i,2)*b(2,n-2)
3912 c(i,n-1) = a(i,1)*b(1,n-1)
3913 $ + a(i,2)*b(2,n-1)
3914 c(i,n) = a(i,1)*b(1,n)
3915 $ + a(i,2)*b(2,n)
3916 enddo
3917 endif
3918
3919 return
3920 end
3921c-----------------------------------------------------------------------
3922 subroutine mxmtest(s,nn,cn,mxmt,name,k,ivb)
3923
3924 real s(nn,2) ! MFLOPS
3925 character*5 cn ! name
3926 character*5 name
3927 external mxmt
3928
3929 include 'SIZE'
3930 parameter (lt=4*lx1*ly1*lz1*lelt)
3931 common /scrns/ a(lt)
3932 common /scruz/ b(lt)
3933 common /scrmg/ c(lt)
3934
3935 integer ll,icalld
3936 save ll,icalld
3937 data ll,icalld /1,0/
3938
3939 if (icalld.eq.0) then ! Initialize matrices:
3940 icalld = icalld + 1
3941 time1 = dnekclock()
3942 call initab(a,b,lt)
3943 time2 = dnekclock()-time1
3944 if (nid.eq.0) write(6,*) 'mxm test init:',lt,time2,name
3945 endif
3946
3947
3948 cn = name
3949
3950c Rectangular matrix tests
3951
3952 nn0 = 1
3953 nn1 = nn
3954 if (ivb.eq.0) then
3955 nn0 = lx1
3956 nn1 = lx1
3957 endif
3958
3959 m = k
3960 do n=nn0,nn1
3961 n1 = n
3962 n2 = n
3963 n3 = n
3964 if (m.eq.1) n1 = n*n
3965 if (m.eq.3) n3 = n*n
3966 if (lt.gt.n1*n3) then
3967 n13 = max(n1,n3)
3968 loop = 250000/(n1*n2*n3) + 500
3969 if (name.eq.'madd ') loop = 200000/(n1*n3) + 5000
3970
3971c-------------------------------------------------------
3972c mem test
3973c-------------------------------------------------------
3974
3975 t0 = dnekclock()
3976 overh = dnekclock()-t0
3977 time1 = dnekclock()
3978 do l=1,loop
3979 if (ll.ge.lt-n1*n3) ll = 1
3980 call mxmt(a(ll),n1,b(ll),n2,c(ll),n3)
3981 ll = ll+n1*n3
3982 enddo
3983 time2 = dnekclock()
3984 time = time2-time1 - overh
3985 iops=loop*n1*n3*(2*n2-1)
3986 if (name.eq.'madd ') iops = loop*n1*n3
3987c write(6,*) loop,time,time2,time1,overh
3988 flops=iops/(1.0e6*time)
3989 s(n,1) = flops
3990c
3991 timel = time/loop
3992 if (nid.eq.0) write(6,199) n,n1,n2,n3,flops,timel,name
3993 199 format(i3,'m',1x,3i6,f10.4,e16.5,3x,a5,' mem')
3994c
3995c-------------------------------------------------------
3996c fast test
3997c-------------------------------------------------------
3998c
3999 call mxmt(a,n1,b,n2,c,n3)
4000 t0 = dnekclock()
4001 overh = dnekclock()-t0
4002 time1 = dnekclock()
4003 do l=1,loop
4004 call mxmt(a,n1,b,n2,c,n3)
4005 enddo
4006 time2 = dnekclock()
4007 time = time2-time1 - overh
4008 iops=loop*n1*n3*(2*n2-1)
4009 if (name.eq.'madd ') iops = loop*n1*n3
4010 flops=iops/(1.0e6*time)
4011 s(n,2) = flops
4012 timel = time/loop
4013c
4014 if (nid.eq.0) write(6,198) n,n1,n2,n3,flops,timel,name
4015 198 format(i3,'f',1x,3i6,f10.4,e16.5,3x,a5,' fast')
4016c
4017 endif
4018 enddo
4019c
4020 return
4021 end
4022c-----------------------------------------------------------------------
4023 subroutine mxm_analyze(s,a,nn,c,nt,ivb)
4024 include 'SIZE'
4025
4026 character*5 c(3,nt)
4027 real s(nn,2,nt,3) ! Measured Mflops, 3 cases
4028 real a(nn,2,nt,3)
4029c ^ ^ ^ |__ N^2xN, NxN, NxN^2
4030c matrix order N __| | |__________which mxm
4031c |
4032c |__cached vs. noncached data
4033
4034
4035 integer itmax(200)
4036
4037 nn0 = 1
4038 nn1 = nn
4039 if (ivb.eq.0) then
4040 nn0 = lx1
4041 nn1 = lx1
4042 endif
4043
4044 do n = nn0,nn1
4045 fmax = 0. ! Peak mflops
4046 do it=1,nt
4047 ai = 0.
4048 di = 0.
4049 do k=1,3
4050 if (s(n,1,it,k).gt.0) then ! Take harmonic means of
4051 ai = ai + 1./s(n,1,it,k) ! case I II and III for
4052 di = di + 1. ! mem test, s(n,1...).
4053 endif
4054 enddo
4055 if (ai.gt.0) ai = di/ai
4056 a(n,1,it,1) = di/ai
4057 if (ai.gt.fmax.and.c(2,it).ne.'madd ') then
4058 fmax = ai
4059 itmax(n) = it
4060 endif
4061 enddo
4062 it = itmax(n)
4063 if (nid.eq.0) write(6,3) n,it,c(2,it),(s(n,1,it,k),k=1,3),fmax
4064 3 format(i3,i2,1x,a5,4f12.0,' Peak harmonic')
4065 enddo
4066 call out_anal(s,a,nn,c,nt,itmax,'Harmonic',1,ivb)
4067c
4068c Case by case
4069c
4070 do k=1,3
4071 do n = nn0,nn1
4072 fmax = 0. ! Peak mflops
4073 do it=1,nt
4074 ai = s(n,1,it,k)
4075 if (ai.gt.fmax.and.c(2,it).ne.'madd ') then
4076 fmax = ai
4077 itmax(n) = it
4078 endif
4079 enddo
4080 enddo
4081 if (k.eq.1) call out_anal(s,a,nn,c,nt,itmax,'Case N2N',k,ivb)
4082 if (k.eq.2) call out_anal(s,a,nn,c,nt,itmax,'Case NxN',k,ivb)
4083 if (k.eq.3) call out_anal(s,a,nn,c,nt,itmax,'Case NN2',k,ivb)
4084 enddo
4085
4086 return
4087 end
4088c-----------------------------------------------------------------------
4089 subroutine out_anal(s,a,nn,c,nt,itmax,name8,k,ivb)
4090 include 'SIZE'
4091
4092 character*5 c(3,nt)
4093 real s(nn,2,nt,3)
4094 real a(nn,2,nt,3)
4095 integer itmax(200)
4096 character*8 name8
4097
4098 if (nid.ne.0) return
4099
4100 nn0 = 1
4101 nn1 = nn
4102 if (ivb.eq.0) then
4103 nn0 = lx1
4104 nn1 = lx1
4105 endif
4106
4107
4108 do n=nn0,nn1
4109 it = itmax(n)
4110 write(6,1) n,s(n,1,it,k),c(2,it),name8
4111 1 format(i4,f14.0,4x,a5,4x,a8,' MxM MFLOPS')
4112 enddo
4113
4114 return
4115 end
4116c-----------------------------------------------------------------------
4117 subroutine mxma(a,n1,b,n2,c,n3)
4118C
4119C Compute C = A*B , for contiguously packed matrices A,B, and C.
4120C
4121C
4122C Compile with -r8 option for 64-bit arithmetic, or convert real
4123c to real*8
4124c
4125 real a(n1,n2),b(n2,n3),c(n1,n3)
4126c
4127 include 'SIZE'
4128 include 'OPCTR'
4129 include 'TOTAL'
4130c
4131#ifdef TIMER
4132 if (isclld.eq.0) then
4133 isclld=1
4134 nrout=nrout+1
4135 myrout=nrout
4136 rname(myrout) = 'mxma '
4137 endif
4138 isbcnt = n1*n3*(2*n2-1)
4139 dct(myrout) = dct(myrout) + (isbcnt)
4140 ncall(myrout) = ncall(myrout) + 1
4141 dcount = dcount + (isbcnt)
4142#endif
4143c
4144c
4145 call mxma2(a,n1,b,n2,c,n3) ! In some cases, this is faster.
4146c
4147 return
4148 end
4149c-----------------------------------------------------------------------
4150 subroutine mxma2(a,n1,b,n2,c,n3)
4151c
4152c unrolled loop routine written by paul fischer
4153c
4154c this version computines C = C + A x B
4155c
4156 real a(n1,n2),b(n2,n3),c(n1,n3)
4157C
4158 integer wdsize
4159 save wdsize
4160 data wdsize/0/
4161c
4162c First call: determine word size for dgemm/sgemm discrimination, below.
4163c
4164 if (wdsize.eq.0) then
4165 one = 1.0
4166 eps = 1.e-12
4167 wdsize = 8
4168 if (one+eps.eq.1.0) wdsize = 4
4169 endif
4170c
4171 if (n2.le.8) then
4172 if (n2.eq.1) then
4173 call mxa1(a,n1,b,n2,c,n3)
4174 elseif (n2.eq.2) then
4175 call mxa2(a,n1,b,n2,c,n3)
4176 elseif (n2.eq.3) then
4177 call mxa3(a,n1,b,n2,c,n3)
4178 elseif (n2.eq.4) then
4179 call mxa4(a,n1,b,n2,c,n3)
4180 elseif (n2.eq.5) then
4181 call mxa5(a,n1,b,n2,c,n3)
4182 elseif (n2.eq.6) then
4183 call mxa6(a,n1,b,n2,c,n3)
4184 elseif (n2.eq.7) then
4185 call mxa7(a,n1,b,n2,c,n3)
4186 else
4187 call mxa8(a,n1,b,n2,c,n3)
4188 endif
4189 elseif (n2.le.16) then
4190 if (n2.eq.9) then
4191 call mxa9(a,n1,b,n2,c,n3)
4192 elseif (n2.eq.10) then
4193 call mxa10(a,n1,b,n2,c,n3)
4194 elseif (n2.eq.11) then
4195 call mxa11(a,n1,b,n2,c,n3)
4196 elseif (n2.eq.12) then
4197 call mxa12(a,n1,b,n2,c,n3)
4198 elseif (n2.eq.13) then
4199 call mxa13(a,n1,b,n2,c,n3)
4200 elseif (n2.eq.14) then
4201 call mxa14(a,n1,b,n2,c,n3)
4202 elseif (n2.eq.15) then
4203 call mxa15(a,n1,b,n2,c,n3)
4204 else
4205 call mxa16(a,n1,b,n2,c,n3)
4206 endif
4207 elseif (n2.le.24) then
4208 if (n2.eq.17) then
4209 call mxa17(a,n1,b,n2,c,n3)
4210 elseif (n2.eq.18) then
4211 call mxa18(a,n1,b,n2,c,n3)
4212 elseif (n2.eq.19) then
4213 call mxa19(a,n1,b,n2,c,n3)
4214 elseif (n2.eq.20) then
4215 call mxa20(a,n1,b,n2,c,n3)
4216 elseif (n2.eq.21) then
4217 call mxa21(a,n1,b,n2,c,n3)
4218 elseif (n2.eq.22) then
4219 call mxa22(a,n1,b,n2,c,n3)
4220 elseif (n2.eq.23) then
4221 call mxa23(a,n1,b,n2,c,n3)
4222 elseif (n2.eq.24) then
4223 call mxa24(a,n1,b,n2,c,n3)
4224 endif
4225 else
4226 call mxm44_0a(a,n1,b,n2,c,n3)
4227 endif
4228c
4229 return
4230 end
4231c-----------------------------------------------------------------------
4232 subroutine mxa1(a,n1,b,n2,c,n3)
4233c
4234 real a(n1,1),b(1,n3),c(n1,n3)
4235c
4236 do j=1,n3
4237 do i=1,n1
4238 c(i,j) = c(i,j)
4239 $ + a(i,1)*b(1,j)
4240 enddo
4241 enddo
4242 return
4243 end
4244c-----------------------------------------------------------------------
4245 subroutine mxa2(a,n1,b,n2,c,n3)
4246c
4247 real a(n1,2),b(2,n3),c(n1,n3)
4248c
4249 do j=1,n3
4250 do i=1,n1
4251 c(i,j) = c(i,j)
4252 $ + a(i,1)*b(1,j)
4253 $ + a(i,2)*b(2,j)
4254 enddo
4255 enddo
4256 return
4257 end
4258c-----------------------------------------------------------------------
4259 subroutine mxa3(a,n1,b,n2,c,n3)
4260c
4261 real a(n1,3),b(3,n3),c(n1,n3)
4262c
4263 do j=1,n3
4264 do i=1,n1
4265 c(i,j) = c(i,j)
4266 $ + a(i,1)*b(1,j)
4267 $ + a(i,2)*b(2,j)
4268 $ + a(i,3)*b(3,j)
4269 enddo
4270 enddo
4271 return
4272 end
4273c-----------------------------------------------------------------------
4274 subroutine mxa4(a,n1,b,n2,c,n3)
4275c
4276 real a(n1,4),b(4,n3),c(n1,n3)
4277c
4278 do j=1,n3
4279 do i=1,n1
4280 c(i,j) = c(i,j)
4281 $ + a(i,1)*b(1,j)
4282 $ + a(i,2)*b(2,j)
4283 $ + a(i,3)*b(3,j)
4284 $ + a(i,4)*b(4,j)
4285 enddo
4286 enddo
4287 return
4288 end
4289c-----------------------------------------------------------------------
4290 subroutine mxa5(a,n1,b,n2,c,n3)
4291c
4292 real a(n1,5),b(5,n3),c(n1,n3)
4293c
4294 do j=1,n3
4295 do i=1,n1
4296 c(i,j) = c(i,j)
4297 $ + a(i,1)*b(1,j)
4298 $ + a(i,2)*b(2,j)
4299 $ + a(i,3)*b(3,j)
4300 $ + a(i,4)*b(4,j)
4301 $ + a(i,5)*b(5,j)
4302 enddo
4303 enddo
4304 return
4305 end
4306c-----------------------------------------------------------------------
4307 subroutine mxa6(a,n1,b,n2,c,n3)
4308c
4309 real a(n1,6),b(6,n3),c(n1,n3)
4310c
4311 do j=1,n3
4312 do i=1,n1
4313 c(i,j) = c(i,j)
4314 $ + a(i,1)*b(1,j)
4315 $ + a(i,2)*b(2,j)
4316 $ + a(i,3)*b(3,j)
4317 $ + a(i,4)*b(4,j)
4318 $ + a(i,5)*b(5,j)
4319 $ + a(i,6)*b(6,j)
4320 enddo
4321 enddo
4322 return
4323 end
4324c-----------------------------------------------------------------------
4325 subroutine mxa7(a,n1,b,n2,c,n3)
4326c
4327 real a(n1,7),b(7,n3),c(n1,n3)
4328c
4329 do j=1,n3
4330 do i=1,n1
4331 c(i,j) = c(i,j)
4332 $ + a(i,1)*b(1,j)
4333 $ + a(i,2)*b(2,j)
4334 $ + a(i,3)*b(3,j)
4335 $ + a(i,4)*b(4,j)
4336 $ + a(i,5)*b(5,j)
4337 $ + a(i,6)*b(6,j)
4338 $ + a(i,7)*b(7,j)
4339 enddo
4340 enddo
4341 return
4342 end
4343c-----------------------------------------------------------------------
4344 subroutine mxa8(a,n1,b,n2,c,n3)
4345c
4346 real a(n1,8),b(8,n3),c(n1,n3)
4347c
4348 do j=1,n3
4349 do i=1,n1
4350 c(i,j) = c(i,j)
4351 $ + a(i,1)*b(1,j)
4352 $ + a(i,2)*b(2,j)
4353 $ + a(i,3)*b(3,j)
4354 $ + a(i,4)*b(4,j)
4355 $ + a(i,5)*b(5,j)
4356 $ + a(i,6)*b(6,j)
4357 $ + a(i,7)*b(7,j)
4358 $ + a(i,8)*b(8,j)
4359 enddo
4360 enddo
4361 return
4362 end
4363c-----------------------------------------------------------------------
4364 subroutine mxa9(a,n1,b,n2,c,n3)
4365c
4366 real a(n1,9),b(9,n3),c(n1,n3)
4367c
4368 do j=1,n3
4369 do i=1,n1
4370 c(i,j) = c(i,j)
4371 $ + a(i,1)*b(1,j)
4372 $ + a(i,2)*b(2,j)
4373 $ + a(i,3)*b(3,j)
4374 $ + a(i,4)*b(4,j)
4375 $ + a(i,5)*b(5,j)
4376 $ + a(i,6)*b(6,j)
4377 $ + a(i,7)*b(7,j)
4378 $ + a(i,8)*b(8,j)
4379 $ + a(i,9)*b(9,j)
4380 enddo
4381 enddo
4382 return
4383 end
4384c-----------------------------------------------------------------------
4385 subroutine mxa10(a,n1,b,n2,c,n3)
4386c
4387 real a(n1,10),b(10,n3),c(n1,n3)
4388c
4389 do j=1,n3
4390 do i=1,n1
4391 c(i,j) = c(i,j)
4392 $ + a(i,1)*b(1,j)
4393 $ + a(i,2)*b(2,j)
4394 $ + a(i,3)*b(3,j)
4395 $ + a(i,4)*b(4,j)
4396 $ + a(i,5)*b(5,j)
4397 $ + a(i,6)*b(6,j)
4398 $ + a(i,7)*b(7,j)
4399 $ + a(i,8)*b(8,j)
4400 $ + a(i,9)*b(9,j)
4401 $ + a(i,10)*b(10,j)
4402 enddo
4403 enddo
4404 return
4405 end
4406c-----------------------------------------------------------------------
4407 subroutine mxa11(a,n1,b,n2,c,n3)
4408c
4409 real a(n1,11),b(11,n3),c(n1,n3)
4410c
4411 do j=1,n3
4412 do i=1,n1
4413 c(i,j) = c(i,j)
4414 $ + a(i,1)*b(1,j)
4415 $ + a(i,2)*b(2,j)
4416 $ + a(i,3)*b(3,j)
4417 $ + a(i,4)*b(4,j)
4418 $ + a(i,5)*b(5,j)
4419 $ + a(i,6)*b(6,j)
4420 $ + a(i,7)*b(7,j)
4421 $ + a(i,8)*b(8,j)
4422 $ + a(i,9)*b(9,j)
4423 $ + a(i,10)*b(10,j)
4424 $ + a(i,11)*b(11,j)
4425 enddo
4426 enddo
4427 return
4428 end
4429c-----------------------------------------------------------------------
4430 subroutine mxa12(a,n1,b,n2,c,n3)
4431c
4432 real a(n1,12),b(12,n3),c(n1,n3)
4433c
4434 do j=1,n3
4435 do i=1,n1
4436 c(i,j) = c(i,j)
4437 $ + a(i,1)*b(1,j)
4438 $ + a(i,2)*b(2,j)
4439 $ + a(i,3)*b(3,j)
4440 $ + a(i,4)*b(4,j)
4441 $ + a(i,5)*b(5,j)
4442 $ + a(i,6)*b(6,j)
4443 $ + a(i,7)*b(7,j)
4444 $ + a(i,8)*b(8,j)
4445 $ + a(i,9)*b(9,j)
4446 $ + a(i,10)*b(10,j)
4447 $ + a(i,11)*b(11,j)
4448 $ + a(i,12)*b(12,j)
4449 enddo
4450 enddo
4451 return
4452 end
4453c-----------------------------------------------------------------------
4454 subroutine mxa13(a,n1,b,n2,c,n3)
4455c
4456 real a(n1,13),b(13,n3),c(n1,n3)
4457c
4458 do j=1,n3
4459 do i=1,n1
4460 c(i,j) = c(i,j)
4461 $ + a(i,1)*b(1,j)
4462 $ + a(i,2)*b(2,j)
4463 $ + a(i,3)*b(3,j)
4464 $ + a(i,4)*b(4,j)
4465 $ + a(i,5)*b(5,j)
4466 $ + a(i,6)*b(6,j)
4467 $ + a(i,7)*b(7,j)
4468 $ + a(i,8)*b(8,j)
4469 $ + a(i,9)*b(9,j)
4470 $ + a(i,10)*b(10,j)
4471 $ + a(i,11)*b(11,j)
4472 $ + a(i,12)*b(12,j)
4473 $ + a(i,13)*b(13,j)
4474 enddo
4475 enddo
4476 return
4477 end
4478c-----------------------------------------------------------------------
4479 subroutine mxa14(a,n1,b,n2,c,n3)
4480c
4481 real a(n1,14),b(14,n3),c(n1,n3)
4482c
4483 do j=1,n3
4484 do i=1,n1
4485 c(i,j) = c(i,j)
4486 $ + a(i,1)*b(1,j)
4487 $ + a(i,2)*b(2,j)
4488 $ + a(i,3)*b(3,j)
4489 $ + a(i,4)*b(4,j)
4490 $ + a(i,5)*b(5,j)
4491 $ + a(i,6)*b(6,j)
4492 $ + a(i,7)*b(7,j)
4493 $ + a(i,8)*b(8,j)
4494 $ + a(i,9)*b(9,j)
4495 $ + a(i,10)*b(10,j)
4496 $ + a(i,11)*b(11,j)
4497 $ + a(i,12)*b(12,j)
4498 $ + a(i,13)*b(13,j)
4499 $ + a(i,14)*b(14,j)
4500 enddo
4501 enddo
4502 return
4503 end
4504c-----------------------------------------------------------------------
4505 subroutine mxa15(a,n1,b,n2,c,n3)
4506c
4507 real a(n1,15),b(15,n3),c(n1,n3)
4508c
4509 do j=1,n3
4510 do i=1,n1
4511 c(i,j) = c(i,j)
4512 $ + a(i,1)*b(1,j)
4513 $ + a(i,2)*b(2,j)
4514 $ + a(i,3)*b(3,j)
4515 $ + a(i,4)*b(4,j)
4516 $ + a(i,5)*b(5,j)
4517 $ + a(i,6)*b(6,j)
4518 $ + a(i,7)*b(7,j)
4519 $ + a(i,8)*b(8,j)
4520 $ + a(i,9)*b(9,j)
4521 $ + a(i,10)*b(10,j)
4522 $ + a(i,11)*b(11,j)
4523 $ + a(i,12)*b(12,j)
4524 $ + a(i,13)*b(13,j)
4525 $ + a(i,14)*b(14,j)
4526 $ + a(i,15)*b(15,j)
4527 enddo
4528 enddo
4529 return
4530 end
4531c-----------------------------------------------------------------------
4532 subroutine mxa16(a,n1,b,n2,c,n3)
4533c
4534 real a(n1,16),b(16,n3),c(n1,n3)
4535c
4536 do j=1,n3
4537 do i=1,n1
4538 c(i,j) = c(i,j)
4539 $ + a(i,1)*b(1,j)
4540 $ + a(i,2)*b(2,j)
4541 $ + a(i,3)*b(3,j)
4542 $ + a(i,4)*b(4,j)
4543 $ + a(i,5)*b(5,j)
4544 $ + a(i,6)*b(6,j)
4545 $ + a(i,7)*b(7,j)
4546 $ + a(i,8)*b(8,j)
4547 $ + a(i,9)*b(9,j)
4548 $ + a(i,10)*b(10,j)
4549 $ + a(i,11)*b(11,j)
4550 $ + a(i,12)*b(12,j)
4551 $ + a(i,13)*b(13,j)
4552 $ + a(i,14)*b(14,j)
4553 $ + a(i,15)*b(15,j)
4554 $ + a(i,16)*b(16,j)
4555 enddo
4556 enddo
4557 return
4558 end
4559c-----------------------------------------------------------------------
4560 subroutine mxa17(a,n1,b,n2,c,n3)
4561c
4562 real a(n1,17),b(17,n3),c(n1,n3)
4563c
4564 do j=1,n3
4565 do i=1,n1
4566 c(i,j) = c(i,j)
4567 $ + a(i,1)*b(1,j)
4568 $ + a(i,2)*b(2,j)
4569 $ + a(i,3)*b(3,j)
4570 $ + a(i,4)*b(4,j)
4571 $ + a(i,5)*b(5,j)
4572 $ + a(i,6)*b(6,j)
4573 $ + a(i,7)*b(7,j)
4574 $ + a(i,8)*b(8,j)
4575 $ + a(i,9)*b(9,j)
4576 $ + a(i,10)*b(10,j)
4577 $ + a(i,11)*b(11,j)
4578 $ + a(i,12)*b(12,j)
4579 $ + a(i,13)*b(13,j)
4580 $ + a(i,14)*b(14,j)
4581 $ + a(i,15)*b(15,j)
4582 $ + a(i,16)*b(16,j)
4583 $ + a(i,17)*b(17,j)
4584 enddo
4585 enddo
4586 return
4587 end
4588c-----------------------------------------------------------------------
4589 subroutine mxa18(a,n1,b,n2,c,n3)
4590c
4591 real a(n1,18),b(18,n3),c(n1,n3)
4592c
4593 do j=1,n3
4594 do i=1,n1
4595 c(i,j) = c(i,j)
4596 $ + a(i,1)*b(1,j)
4597 $ + a(i,2)*b(2,j)
4598 $ + a(i,3)*b(3,j)
4599 $ + a(i,4)*b(4,j)
4600 $ + a(i,5)*b(5,j)
4601 $ + a(i,6)*b(6,j)
4602 $ + a(i,7)*b(7,j)
4603 $ + a(i,8)*b(8,j)
4604 $ + a(i,9)*b(9,j)
4605 $ + a(i,10)*b(10,j)
4606 $ + a(i,11)*b(11,j)
4607 $ + a(i,12)*b(12,j)
4608 $ + a(i,13)*b(13,j)
4609 $ + a(i,14)*b(14,j)
4610 $ + a(i,15)*b(15,j)
4611 $ + a(i,16)*b(16,j)
4612 $ + a(i,17)*b(17,j)
4613 $ + a(i,18)*b(18,j)
4614 enddo
4615 enddo
4616 return
4617 end
4618c-----------------------------------------------------------------------
4619 subroutine mxa19(a,n1,b,n2,c,n3)
4620c
4621 real a(n1,19),b(19,n3),c(n1,n3)
4622c
4623 do j=1,n3
4624 do i=1,n1
4625 c(i,j) = c(i,j)
4626 $ + a(i,1)*b(1,j)
4627 $ + a(i,2)*b(2,j)
4628 $ + a(i,3)*b(3,j)
4629 $ + a(i,4)*b(4,j)
4630 $ + a(i,5)*b(5,j)
4631 $ + a(i,6)*b(6,j)
4632 $ + a(i,7)*b(7,j)
4633 $ + a(i,8)*b(8,j)
4634 $ + a(i,9)*b(9,j)
4635 $ + a(i,10)*b(10,j)
4636 $ + a(i,11)*b(11,j)
4637 $ + a(i,12)*b(12,j)
4638 $ + a(i,13)*b(13,j)
4639 $ + a(i,14)*b(14,j)
4640 $ + a(i,15)*b(15,j)
4641 $ + a(i,16)*b(16,j)
4642 $ + a(i,17)*b(17,j)
4643 $ + a(i,18)*b(18,j)
4644 $ + a(i,19)*b(19,j)
4645 enddo
4646 enddo
4647 return
4648 end
4649c-----------------------------------------------------------------------
4650 subroutine mxa20(a,n1,b,n2,c,n3)
4651c
4652 real a(n1,20),b(20,n3),c(n1,n3)
4653c
4654 do j=1,n3
4655 do i=1,n1
4656 c(i,j) = c(i,j)
4657 $ + a(i,1)*b(1,j)
4658 $ + a(i,2)*b(2,j)
4659 $ + a(i,3)*b(3,j)
4660 $ + a(i,4)*b(4,j)
4661 $ + a(i,5)*b(5,j)
4662 $ + a(i,6)*b(6,j)
4663 $ + a(i,7)*b(7,j)
4664 $ + a(i,8)*b(8,j)
4665 $ + a(i,9)*b(9,j)
4666 $ + a(i,10)*b(10,j)
4667 $ + a(i,11)*b(11,j)
4668 $ + a(i,12)*b(12,j)
4669 $ + a(i,13)*b(13,j)
4670 $ + a(i,14)*b(14,j)
4671 $ + a(i,15)*b(15,j)
4672 $ + a(i,16)*b(16,j)
4673 $ + a(i,17)*b(17,j)
4674 $ + a(i,18)*b(18,j)
4675 $ + a(i,19)*b(19,j)
4676 $ + a(i,20)*b(20,j)
4677 enddo
4678 enddo
4679 return
4680 end
4681c-----------------------------------------------------------------------
4682 subroutine mxa21(a,n1,b,n2,c,n3)
4683c
4684 real a(n1,21),b(21,n3),c(n1,n3)
4685c
4686 do j=1,n3
4687 do i=1,n1
4688 c(i,j) = c(i,j)
4689 $ + a(i,1)*b(1,j)
4690 $ + a(i,2)*b(2,j)
4691 $ + a(i,3)*b(3,j)
4692 $ + a(i,4)*b(4,j)
4693 $ + a(i,5)*b(5,j)
4694 $ + a(i,6)*b(6,j)
4695 $ + a(i,7)*b(7,j)
4696 $ + a(i,8)*b(8,j)
4697 $ + a(i,9)*b(9,j)
4698 $ + a(i,10)*b(10,j)
4699 $ + a(i,11)*b(11,j)
4700 $ + a(i,12)*b(12,j)
4701 $ + a(i,13)*b(13,j)
4702 $ + a(i,14)*b(14,j)
4703 $ + a(i,15)*b(15,j)
4704 $ + a(i,16)*b(16,j)
4705 $ + a(i,17)*b(17,j)
4706 $ + a(i,18)*b(18,j)
4707 $ + a(i,19)*b(19,j)
4708 $ + a(i,20)*b(20,j)
4709 $ + a(i,21)*b(21,j)
4710 enddo
4711 enddo
4712 return
4713 end
4714c-----------------------------------------------------------------------
4715 subroutine mxa22(a,n1,b,n2,c,n3)
4716c
4717 real a(n1,22),b(22,n3),c(n1,n3)
4718c
4719 do j=1,n3
4720 do i=1,n1
4721 c(i,j) = c(i,j)
4722 $ + a(i,1)*b(1,j)
4723 $ + a(i,2)*b(2,j)
4724 $ + a(i,3)*b(3,j)
4725 $ + a(i,4)*b(4,j)
4726 $ + a(i,5)*b(5,j)
4727 $ + a(i,6)*b(6,j)
4728 $ + a(i,7)*b(7,j)
4729 $ + a(i,8)*b(8,j)
4730 $ + a(i,9)*b(9,j)
4731 $ + a(i,10)*b(10,j)
4732 $ + a(i,11)*b(11,j)
4733 $ + a(i,12)*b(12,j)
4734 $ + a(i,13)*b(13,j)
4735 $ + a(i,14)*b(14,j)
4736 $ + a(i,15)*b(15,j)
4737 $ + a(i,16)*b(16,j)
4738 $ + a(i,17)*b(17,j)
4739 $ + a(i,18)*b(18,j)
4740 $ + a(i,19)*b(19,j)
4741 $ + a(i,20)*b(20,j)
4742 $ + a(i,21)*b(21,j)
4743 $ + a(i,22)*b(22,j)
4744 enddo
4745 enddo
4746 return
4747 end
4748c-----------------------------------------------------------------------
4749 subroutine mxa23(a,n1,b,n2,c,n3)
4750c
4751 real a(n1,23),b(23,n3),c(n1,n3)
4752c
4753 do j=1,n3
4754 do i=1,n1
4755 c(i,j) = c(i,j)
4756 $ + a(i,1)*b(1,j)
4757 $ + a(i,2)*b(2,j)
4758 $ + a(i,3)*b(3,j)
4759 $ + a(i,4)*b(4,j)
4760 $ + a(i,5)*b(5,j)
4761 $ + a(i,6)*b(6,j)
4762 $ + a(i,7)*b(7,j)
4763 $ + a(i,8)*b(8,j)
4764 $ + a(i,9)*b(9,j)
4765 $ + a(i,10)*b(10,j)
4766 $ + a(i,11)*b(11,j)
4767 $ + a(i,12)*b(12,j)
4768 $ + a(i,13)*b(13,j)
4769 $ + a(i,14)*b(14,j)
4770 $ + a(i,15)*b(15,j)
4771 $ + a(i,16)*b(16,j)
4772 $ + a(i,17)*b(17,j)
4773 $ + a(i,18)*b(18,j)
4774 $ + a(i,19)*b(19,j)
4775 $ + a(i,20)*b(20,j)
4776 $ + a(i,21)*b(21,j)
4777 $ + a(i,22)*b(22,j)
4778 $ + a(i,23)*b(23,j)
4779 enddo
4780 enddo
4781 return
4782 end
4783c-----------------------------------------------------------------------
4784 subroutine mxa24(a,n1,b,n2,c,n3)
4785c
4786 real a(n1,24),b(24,n3),c(n1,n3)
4787c
4788 do j=1,n3
4789 do i=1,n1
4790 c(i,j) = c(i,j)
4791 $ + a(i,1)*b(1,j)
4792 $ + a(i,2)*b(2,j)
4793 $ + a(i,3)*b(3,j)
4794 $ + a(i,4)*b(4,j)
4795 $ + a(i,5)*b(5,j)
4796 $ + a(i,6)*b(6,j)
4797 $ + a(i,7)*b(7,j)
4798 $ + a(i,8)*b(8,j)
4799 $ + a(i,9)*b(9,j)
4800 $ + a(i,10)*b(10,j)
4801 $ + a(i,11)*b(11,j)
4802 $ + a(i,12)*b(12,j)
4803 $ + a(i,13)*b(13,j)
4804 $ + a(i,14)*b(14,j)
4805 $ + a(i,15)*b(15,j)
4806 $ + a(i,16)*b(16,j)
4807 $ + a(i,17)*b(17,j)
4808 $ + a(i,18)*b(18,j)
4809 $ + a(i,19)*b(19,j)
4810 $ + a(i,20)*b(20,j)
4811 $ + a(i,21)*b(21,j)
4812 $ + a(i,22)*b(22,j)
4813 $ + a(i,23)*b(23,j)
4814 $ + a(i,24)*b(24,j)
4815 enddo
4816 enddo
4817 return
4818 end
4819c-----------------------------------------------------------------------
4820 subroutine mxm44_0a(a, m, b, k, c, n)
4821c
4822c routine written by Bruce Curtiss
4823c
4824c subroutine matmul44(m, n, k, a, lda, b, ldb, c, ldc)
4825c real a(lda,k), b(ldb,n), c(ldc,n)
4826c
4827 real a(m,k), b(k,n), c(m,n)
4828 real s11, s12, s13, s14, s21, s22, s23, s24
4829 real s31, s32, s33, s34, s41, s42, s43, s44
4830c
4831c matrix multiply with a 4x4 pencil
4832c
4833
4834 mresid = mod(m,4)
4835 nresid = mod(n,4)
4836 m1 = m - mresid + 1
4837 n1 = n - nresid + 1
4838
4839 do i=1,m-mresid,4
4840 do j=1,n-nresid,4
4841 s11 = c(i,j)
4842 s12 = c(i,j+1)
4843 s13 = c(i,j+2)
4844 s14 = c(i,j+3)
4845
4846 s21 = c(i+1,j)
4847 s31 = c(i+2,j)
4848 s41 = c(i+3,j)
4849
4850 s22 = c(i+1,j+1)
4851 s32 = c(i+2,j+1)
4852 s42 = c(i+3,j+1)
4853
4854 s23 = c(i+1,j+2)
4855 s33 = c(i+2,j+2)
4856 s43 = c(i+3,j+2)
4857
4858 s24 = c(i+1,j+3)
4859 s34 = c(i+2,j+3)
4860 s44 = c(i+3,j+3)
4861c
4862 do l=1,k
4863 s11 = s11 + a(i,l)*b(l,j)
4864 s12 = s12 + a(i,l)*b(l,j+1)
4865 s13 = s13 + a(i,l)*b(l,j+2)
4866 s14 = s14 + a(i,l)*b(l,j+3)
4867
4868 s21 = s21 + a(i+1,l)*b(l,j)
4869 s22 = s22 + a(i+1,l)*b(l,j+1)
4870 s23 = s23 + a(i+1,l)*b(l,j+2)
4871 s24 = s24 + a(i+1,l)*b(l,j+3)
4872
4873 s31 = s31 + a(i+2,l)*b(l,j)
4874 s32 = s32 + a(i+2,l)*b(l,j+1)
4875 s33 = s33 + a(i+2,l)*b(l,j+2)
4876 s34 = s34 + a(i+2,l)*b(l,j+3)
4877
4878 s41 = s41 + a(i+3,l)*b(l,j)
4879 s42 = s42 + a(i+3,l)*b(l,j+1)
4880 s43 = s43 + a(i+3,l)*b(l,j+2)
4881 s44 = s44 + a(i+3,l)*b(l,j+3)
4882 enddo
4883 c(i,j) = s11
4884 c(i,j+1) = s12
4885 c(i,j+2) = s13
4886 c(i,j+3) = s14
4887
4888 c(i+1,j) = s21
4889 c(i+2,j) = s31
4890 c(i+3,j) = s41
4891
4892 c(i+1,j+1) = s22
4893 c(i+2,j+1) = s32
4894 c(i+3,j+1) = s42
4895
4896 c(i+1,j+2) = s23
4897 c(i+2,j+2) = s33
4898 c(i+3,j+2) = s43
4899
4900 c(i+1,j+3) = s24
4901 c(i+2,j+3) = s34
4902 c(i+3,j+3) = s44
4903 enddo
4904c Residual when n is not multiple of 4
4905 if (nresid .ne. 0) then
4906 if (nresid .eq. 1) then
4907 s11 = c(i ,n)
4908 s21 = c(i+1,n)
4909 s31 = c(i+2,n)
4910 s41 = c(i+3,n)
4911 do l=1,k
4912 s11 = s11 + a(i,l)*b(l,n)
4913 s21 = s21 + a(i+1,l)*b(l,n)
4914 s31 = s31 + a(i+2,l)*b(l,n)
4915 s41 = s41 + a(i+3,l)*b(l,n)
4916 enddo
4917 c(i,n) = s11
4918 c(i+1,n) = s21
4919 c(i+2,n) = s31
4920 c(i+3,n) = s41
4921 elseif (nresid .eq. 2) then
4922 s11 = c(i ,j )
4923 s12 = c(i ,j+1)
4924 s21 = c(i+1,j )
4925 s31 = c(i+2,j )
4926 s41 = c(i+3,j )
4927 s22 = c(i+1,j+1)
4928 s32 = c(i+2,j+1)
4929 s42 = c(i+3,j+1)
4930 do l=1,k
4931 s11 = s11 + a(i,l)*b(l,j)
4932 s12 = s12 + a(i,l)*b(l,j+1)
4933
4934 s21 = s21 + a(i+1,l)*b(l,j)
4935 s22 = s22 + a(i+1,l)*b(l,j+1)
4936
4937 s31 = s31 + a(i+2,l)*b(l,j)
4938 s32 = s32 + a(i+2,l)*b(l,j+1)
4939
4940 s41 = s41 + a(i+3,l)*b(l,j)
4941 s42 = s42 + a(i+3,l)*b(l,j+1)
4942 enddo
4943 c(i ,j ) = s11
4944 c(i ,j+1) = s12
4945 c(i+1,j ) = s21
4946 c(i+2,j ) = s31
4947 c(i+3,j ) = s41
4948 c(i+1,j+1) = s22
4949 c(i+2,j+1) = s32
4950 c(i+3,j+1) = s42
4951 else
4952 s11 = c(i,j)
4953 s21 = c(i+1,j)
4954 s31 = c(i+2,j)
4955 s41 = c(i+3,j)
4956 s12 = c(i,j+1)
4957 s22 = c(i+1,j+1)
4958 s32 = c(i+2,j+1)
4959 s42 = c(i+3,j+1)
4960 s13 = c(i,j+2)
4961 s23 = c(i+1,j+2)
4962 s33 = c(i+2,j+2)
4963 s43 = c(i+3,j+2)
4964 do l=1,k
4965 s11 = s11 + a(i,l)*b(l,j)
4966 s12 = s12 + a(i,l)*b(l,j+1)
4967 s13 = s13 + a(i,l)*b(l,j+2)
4968
4969 s21 = s21 + a(i+1,l)*b(l,j)
4970 s22 = s22 + a(i+1,l)*b(l,j+1)
4971 s23 = s23 + a(i+1,l)*b(l,j+2)
4972
4973 s31 = s31 + a(i+2,l)*b(l,j)
4974 s32 = s32 + a(i+2,l)*b(l,j+1)
4975 s33 = s33 + a(i+2,l)*b(l,j+2)
4976
4977 s41 = s41 + a(i+3,l)*b(l,j)
4978 s42 = s42 + a(i+3,l)*b(l,j+1)
4979 s43 = s43 + a(i+3,l)*b(l,j+2)
4980 enddo
4981 c(i,j) = s11
4982 c(i+1,j) = s21
4983 c(i+2,j) = s31
4984 c(i+3,j) = s41
4985 c(i,j+1) = s12
4986 c(i+1,j+1) = s22
4987 c(i+2,j+1) = s32
4988 c(i+3,j+1) = s42
4989 c(i,j+2) = s13
4990 c(i+1,j+2) = s23
4991 c(i+2,j+2) = s33
4992 c(i+3,j+2) = s43
4993 endif
4994 endif
4995 enddo
4996
4997c Residual when m is not multiple of 4
4998 if (mresid .eq. 0) then
4999 return
5000 elseif (mresid .eq. 1) then
5001 do j=1,n-nresid,4
5002 s11 = c(m,j)
5003 s12 = c(m,j+1)
5004 s13 = c(m,j+2)
5005 s14 = c(m,j+3)
5006 do l=1,k
5007 s11 = s11 + a(m,l)*b(l,j)
5008 s12 = s12 + a(m,l)*b(l,j+1)
5009 s13 = s13 + a(m,l)*b(l,j+2)
5010 s14 = s14 + a(m,l)*b(l,j+3)
5011 enddo
5012 c(m,j) = s11
5013 c(m,j+1) = s12
5014 c(m,j+2) = s13
5015 c(m,j+3) = s14
5016 enddo
5017c mresid is 1, check nresid
5018 if (nresid .eq. 0) then
5019 return
5020 elseif (nresid .eq. 1) then
5021 s11 = c(m,n)
5022 do l=1,k
5023 s11 = s11 + a(m,l)*b(l,n)
5024 enddo
5025 c(m,n) = s11
5026 return
5027 elseif (nresid .eq. 2) then
5028 s11 = c(m,n-1)
5029 s12 = c(m,n)
5030 do l=1,k
5031 s11 = s11 + a(m,l)*b(l,n-1)
5032 s12 = s12 + a(m,l)*b(l,n)
5033 enddo
5034 c(m,n-1) = s11
5035 c(m,n) = s12
5036 return
5037 else
5038 s11 = c(m,n-2)
5039 s12 = c(m,n-1)
5040 s13 = c(m,n)
5041 do l=1,k
5042 s11 = s11 + a(m,l)*b(l,n-2)
5043 s12 = s12 + a(m,l)*b(l,n-1)
5044 s13 = s13 + a(m,l)*b(l,n)
5045 enddo
5046 c(m,n-2) = s11
5047 c(m,n-1) = s12
5048 c(m,n) = s13
5049 return
5050 endif
5051 elseif (mresid .eq. 2) then
5052 do j=1,n-nresid,4
5053 s11 = c(m-1,j)
5054 s12 = c(m-1,j+1)
5055 s13 = c(m-1,j+2)
5056 s14 = c(m-1,j+3)
5057 s21 = c(m,j)
5058 s22 = c(m,j+1)
5059 s23 = c(m,j+2)
5060 s24 = c(m,j+3)
5061 do l=1,k
5062 s11 = s11 + a(m-1,l)*b(l,j)
5063 s12 = s12 + a(m-1,l)*b(l,j+1)
5064 s13 = s13 + a(m-1,l)*b(l,j+2)
5065 s14 = s14 + a(m-1,l)*b(l,j+3)
5066
5067 s21 = s21 + a(m,l)*b(l,j)
5068 s22 = s22 + a(m,l)*b(l,j+1)
5069 s23 = s23 + a(m,l)*b(l,j+2)
5070 s24 = s24 + a(m,l)*b(l,j+3)
5071 enddo
5072 c(m-1,j) = s11
5073 c(m-1,j+1) = s12
5074 c(m-1,j+2) = s13
5075 c(m-1,j+3) = s14
5076 c(m,j) = s21
5077 c(m,j+1) = s22
5078 c(m,j+2) = s23
5079 c(m,j+3) = s24
5080 enddo
5081c mresid is 2, check nresid
5082 if (nresid .eq. 0) then
5083 return
5084 elseif (nresid .eq. 1) then
5085 s11 = c(m-1,n)
5086 s21 = c(m,n)
5087 do l=1,k
5088 s11 = s11 + a(m-1,l)*b(l,n)
5089 s21 = s21 + a(m,l)*b(l,n)
5090 enddo
5091 c(m-1,n) = s11
5092 c(m,n) = s21
5093 return
5094 elseif (nresid .eq. 2) then
5095 s11 = c(m-1,n-1)
5096 s12 = c(m-1,n)
5097 s21 = c(m,n-1)
5098 s22 = c(m,n)
5099 do l=1,k
5100 s11 = s11 + a(m-1,l)*b(l,n-1)
5101 s12 = s12 + a(m-1,l)*b(l,n)
5102 s21 = s21 + a(m,l)*b(l,n-1)
5103 s22 = s22 + a(m,l)*b(l,n)
5104 enddo
5105 c(m-1,n-1) = s11
5106 c(m-1,n) = s12
5107 c(m,n-1) = s21
5108 c(m,n) = s22
5109 return
5110 else
5111 s11 = c(m-1,n-2)
5112 s12 = c(m-1,n-1)
5113 s13 = c(m-1,n)
5114 s21 = c(m,n-2)
5115 s22 = c(m,n-1)
5116 s23 = c(m,n)
5117 do l=1,k
5118 s11 = s11 + a(m-1,l)*b(l,n-2)
5119 s12 = s12 + a(m-1,l)*b(l,n-1)
5120 s13 = s13 + a(m-1,l)*b(l,n)
5121 s21 = s21 + a(m,l)*b(l,n-2)
5122 s22 = s22 + a(m,l)*b(l,n-1)
5123 s23 = s23 + a(m,l)*b(l,n)
5124 enddo
5125 c(m-1,n-2) = s11
5126 c(m-1,n-1) = s12
5127 c(m-1,n) = s13
5128 c(m,n-2) = s21
5129 c(m,n-1) = s22
5130 c(m,n) = s23
5131 return
5132 endif
5133 else
5134c mresid is 3
5135 do j=1,n-nresid,4
5136 s11 = c(m-2,j)
5137 s12 = c(m-2,j+1)
5138 s13 = c(m-2,j+2)
5139 s14 = c(m-2,j+3)
5140
5141 s21 = c(m-1,j)
5142 s22 = c(m-1,j+1)
5143 s23 = c(m-1,j+2)
5144 s24 = c(m-1,j+3)
5145
5146 s31 = c(m,j)
5147 s32 = c(m,j+1)
5148 s33 = c(m,j+2)
5149 s34 = c(m,j+3)
5150
5151 do l=1,k
5152 s11 = s11 + a(m-2,l)*b(l,j)
5153 s12 = s12 + a(m-2,l)*b(l,j+1)
5154 s13 = s13 + a(m-2,l)*b(l,j+2)
5155 s14 = s14 + a(m-2,l)*b(l,j+3)
5156
5157 s21 = s21 + a(m-1,l)*b(l,j)
5158 s22 = s22 + a(m-1,l)*b(l,j+1)
5159 s23 = s23 + a(m-1,l)*b(l,j+2)
5160 s24 = s24 + a(m-1,l)*b(l,j+3)
5161
5162 s31 = s31 + a(m,l)*b(l,j)
5163 s32 = s32 + a(m,l)*b(l,j+1)
5164 s33 = s33 + a(m,l)*b(l,j+2)
5165 s34 = s34 + a(m,l)*b(l,j+3)
5166 enddo
5167 c(m-2,j) = s11
5168 c(m-2,j+1) = s12
5169 c(m-2,j+2) = s13
5170 c(m-2,j+3) = s14
5171
5172 c(m-1,j) = s21
5173 c(m-1,j+1) = s22
5174 c(m-1,j+2) = s23
5175 c(m-1,j+3) = s24
5176
5177 c(m,j) = s31
5178 c(m,j+1) = s32
5179 c(m,j+2) = s33
5180 c(m,j+3) = s34
5181 enddo
5182c mresid is 3, check nresid
5183 if (nresid .eq. 0) then
5184 return
5185 elseif (nresid .eq. 1) then
5186 s11 = c(m-2,n)
5187 s21 = c(m-1,n)
5188 s31 = c(m,n)
5189 do l=1,k
5190 s11 = s11 + a(m-2,l)*b(l,n)
5191 s21 = s21 + a(m-1,l)*b(l,n)
5192 s31 = s31 + a(m,l)*b(l,n)
5193 enddo
5194 c(m-2,n) = s11
5195 c(m-1,n) = s21
5196 c(m,n) = s31
5197 return
5198 elseif (nresid .eq. 2) then
5199 s11 = c(m-2,n-1)
5200 s12 = c(m-2,n)
5201 s21 = c(m-1,n-1)
5202 s22 = c(m-1,n)
5203 s31 = c(m,n-1)
5204 s32 = c(m,n)
5205 do l=1,k
5206 s11 = s11 + a(m-2,l)*b(l,n-1)
5207 s12 = s12 + a(m-2,l)*b(l,n)
5208 s21 = s21 + a(m-1,l)*b(l,n-1)
5209 s22 = s22 + a(m-1,l)*b(l,n)
5210 s31 = s31 + a(m,l)*b(l,n-1)
5211 s32 = s32 + a(m,l)*b(l,n)
5212 enddo
5213 c(m-2,n-1) = s11
5214 c(m-2,n) = s12
5215 c(m-1,n-1) = s21
5216 c(m-1,n) = s22
5217 c(m,n-1) = s31
5218 c(m,n) = s32
5219 return
5220 else
5221 s11 = c(m-2,n-2)
5222 s12 = c(m-2,n-1)
5223 s13 = c(m-2,n)
5224 s21 = c(m-1,n-2)
5225 s22 = c(m-1,n-1)
5226 s23 = c(m-1,n)
5227 s31 = c(m,n-2)
5228 s32 = c(m,n-1)
5229 s33 = c(m,n)
5230 do l=1,k
5231 s11 = s11 + a(m-2,l)*b(l,n-2)
5232 s12 = s12 + a(m-2,l)*b(l,n-1)
5233 s13 = s13 + a(m-2,l)*b(l,n)
5234 s21 = s21 + a(m-1,l)*b(l,n-2)
5235 s22 = s22 + a(m-1,l)*b(l,n-1)
5236 s23 = s23 + a(m-1,l)*b(l,n)
5237 s31 = s31 + a(m,l)*b(l,n-2)
5238 s32 = s32 + a(m,l)*b(l,n-1)
5239 s33 = s33 + a(m,l)*b(l,n)
5240 enddo
5241 c(m-2,n-2) = s11
5242 c(m-2,n-1) = s12
5243 c(m-2,n) = s13
5244 c(m-1,n-2) = s21
5245 c(m-1,n-1) = s22
5246 c(m-1,n) = s23
5247 c(m,n-2) = s31
5248 c(m,n-1) = s32
5249 c(m,n) = s33
5250 return
5251 endif
5252 endif
5253
5254 return
5255 end
5256c-----------------------------------------------------------------------
5257 subroutine mxm44_2a(a, m, b, k, c, n)
5258 real a(m,2), b(2,n), c(m,n)
5259
5260 nresid = iand(n,3)
5261 n1 = n - nresid + 1
5262
5263 do j=1,n-nresid,4
5264 do i=1,m
5265 c(i,j ) = c(i,j)
5266 $ + a(i,1)*b(1,j)
5267 $ + a(i,2)*b(2,j)
5268 c(i,j+1) = c(i,j+1)
5269 $ + a(i,1)*b(1,j+1)
5270 $ + a(i,2)*b(2,j+1)
5271 c(i,j+2) = c(i,j+2)
5272 $ + a(i,1)*b(1,j+2)
5273 $ + a(i,2)*b(2,j+2)
5274 c(i,j+3) = c(i,j+3)
5275 $ + a(i,1)*b(1,j+3)
5276 $ + a(i,2)*b(2,j+3)
5277 enddo
5278 enddo
5279 if (nresid .eq. 0) then
5280 return
5281 elseif (nresid .eq. 1) then
5282 do i=1,m
5283 c(i,n) = c(i,n)
5284 $ + a(i,1)*b(1,n)
5285 $ + a(i,2)*b(2,n)
5286 enddo
5287 elseif (nresid .eq. 2) then
5288 do i=1,m
5289 c(i,n-1) = c(i,n-1)
5290 $ + a(i,1)*b(1,n-1)
5291 $ + a(i,2)*b(2,n-1)
5292 c(i,n ) = c(i,n )
5293 $ + a(i,1)*b(1,n )
5294 $ + a(i,2)*b(2,n )
5295 enddo
5296 else
5297 do i=1,m
5298 c(i,n-2) = c(i,n-2)
5299 $ + a(i,1)*b(1,n-2)
5300 $ + a(i,2)*b(2,n-2)
5301 c(i,n-1) = c(i,n-1)
5302 $ + a(i,1)*b(1,n-1)
5303 $ + a(i,2)*b(2,n-1)
5304 c(i,n ) = c(i,n )
5305 $ + a(i,1)*b(1,n )
5306 $ + a(i,2)*b(2,n )
5307 enddo
5308 endif
5309
5310 return
5311 end
5312c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.