KTH framework for Nek5000 toolboxes; testing version  0.0.1
mxm_std.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine mxmf2(A,N1,B,N2,C,N3)
3 c
4 c unrolled loop version
5 c
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
65 c
66  return
67  end
68 c-----------------------------------------------------------------------
69  subroutine mxf1(a,n1,b,n2,c,n3)
70 c
71  real a(n1,1),b(1,n3),c(n1,n3)
72 c
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
80 c-----------------------------------------------------------------------
81  subroutine mxf2(a,n1,b,n2,c,n3)
82 c
83  real a(n1,2),b(2,n3),c(n1,n3)
84 c
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
93 c-----------------------------------------------------------------------
94  subroutine mxf3(a,n1,b,n2,c,n3)
95 c
96  real a(n1,3),b(3,n3),c(n1,n3)
97 c
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
107 c-----------------------------------------------------------------------
108  subroutine mxf4(a,n1,b,n2,c,n3)
109 c
110  real a(n1,4),b(4,n3),c(n1,n3)
111 c
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
122 c-----------------------------------------------------------------------
123  subroutine mxf5(a,n1,b,n2,c,n3)
124 c
125  real a(n1,5),b(5,n3),c(n1,n3)
126 c
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
138 c-----------------------------------------------------------------------
139  subroutine mxf6(a,n1,b,n2,c,n3)
140 c
141  real a(n1,6),b(6,n3),c(n1,n3)
142 c
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
155 c-----------------------------------------------------------------------
156  subroutine mxf7(a,n1,b,n2,c,n3)
157 c
158  real a(n1,7),b(7,n3),c(n1,n3)
159 c
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
173 c-----------------------------------------------------------------------
174  subroutine mxf8(a,n1,b,n2,c,n3)
175 c
176  real a(n1,8),b(8,n3),c(n1,n3)
177 c
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
192 c-----------------------------------------------------------------------
193  subroutine mxf9(a,n1,b,n2,c,n3)
194 c
195  real a(n1,9),b(9,n3),c(n1,n3)
196 c
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
212 c-----------------------------------------------------------------------
213  subroutine mxf10(a,n1,b,n2,c,n3)
214 c
215  real a(n1,10),b(10,n3),c(n1,n3)
216 c
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
233 c-----------------------------------------------------------------------
234  subroutine mxf11(a,n1,b,n2,c,n3)
235 c
236  real a(n1,11),b(11,n3),c(n1,n3)
237 c
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
255 c-----------------------------------------------------------------------
256  subroutine mxf12(a,n1,b,n2,c,n3)
257 c
258  real a(n1,12),b(12,n3),c(n1,n3)
259 c
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
278 c-----------------------------------------------------------------------
279  subroutine mxf13(a,n1,b,n2,c,n3)
280 c
281  real a(n1,13),b(13,n3),c(n1,n3)
282 c
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
302 c-----------------------------------------------------------------------
303  subroutine mxf14(a,n1,b,n2,c,n3)
304 c
305  real a(n1,14),b(14,n3),c(n1,n3)
306 c
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
327 c-----------------------------------------------------------------------
328  subroutine mxf15(a,n1,b,n2,c,n3)
329 c
330  real a(n1,15),b(15,n3),c(n1,n3)
331 c
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
353 c-----------------------------------------------------------------------
354  subroutine mxf16(a,n1,b,n2,c,n3)
355 c
356  real a(n1,16),b(16,n3),c(n1,n3)
357 c
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
380 c-----------------------------------------------------------------------
381  subroutine mxf17(a,n1,b,n2,c,n3)
382 c
383  real a(n1,17),b(17,n3),c(n1,n3)
384 c
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
408 c-----------------------------------------------------------------------
409  subroutine mxf18(a,n1,b,n2,c,n3)
410 c
411  real a(n1,18),b(18,n3),c(n1,n3)
412 c
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
437 c-----------------------------------------------------------------------
438  subroutine mxf19(a,n1,b,n2,c,n3)
439 c
440  real a(n1,19),b(19,n3),c(n1,n3)
441 c
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
467 c-----------------------------------------------------------------------
468  subroutine mxf20(a,n1,b,n2,c,n3)
469 c
470  real a(n1,20),b(20,n3),c(n1,n3)
471 c
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
498 c-----------------------------------------------------------------------
499  subroutine mxf21(a,n1,b,n2,c,n3)
500 c
501  real a(n1,21),b(21,n3),c(n1,n3)
502 c
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
530 c-----------------------------------------------------------------------
531  subroutine mxf22(a,n1,b,n2,c,n3)
532 c
533  real a(n1,22),b(22,n3),c(n1,n3)
534 c
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
563 c-----------------------------------------------------------------------
564  subroutine mxf23(a,n1,b,n2,c,n3)
565 c
566  real a(n1,23),b(23,n3),c(n1,n3)
567 c
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
597 c-----------------------------------------------------------------------
598  subroutine mxf24(a,n1,b,n2,c,n3)
599 c
600  real a(n1,24),b(24,n3),c(n1,n3)
601 c
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
632 c-----------------------------------------------------------------------
633  subroutine mxm44_0(a, m, b, k, c, n)
634 c
635 c matrix multiply with a 4x4 pencil
636 c
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
1061 c-----------------------------------------------------------------------
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
1107 c-----------------------------------------------------------------------
1108  subroutine mxm_test_all(nid,ivb)
1109 c
1110 c Collect matrix-matrix product statistics
1111 c
1112  external mxms,mxmur2,mxmur3,mxmd,mxmfb,mxmf3,mxmu4
1113  external madd,mxm,mxm44
1114 c
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
1143 c-----------------------------------------------------------------------
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
1158 c-----------------------------------------------------------------------
1159  subroutine mxms(a,n1,b,n2,c,n3)
1160 C----------------------------------------------------------------------
1161 C
1162 C Matrix-vector product routine.
1163 C NOTE: Use assembly coded routine if available.
1164 C
1165 C---------------------------------------------------------------------
1166  REAL A(N1,N2),B(N2,N3),C(N1,N3)
1167 C
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
1180 c-----------------------------------------------------------------------
1181  subroutine mxmu4(a,n1,b,n2,c,n3)
1182 C----------------------------------------------------------------------
1183 C
1184 C Matrix-vector product routine.
1185 C NOTE: Use assembly coded routine if available.
1186 C
1187 C---------------------------------------------------------------------
1188  REAL A(N1,N2),B(N2,N3),C(N1,N3)
1189 C
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
1210 c-----------------------------------------------------------------------
1211  subroutine madd (a,n1,b,n2,c,n3)
1212 c
1213  real a(n1,n2),b(n2,n3),c(n1,n3)
1214 c
1215  do j=1,n3
1216  do i=1,n1
1217  c(i,j) = a(i,j)+b(i,j)
1218  enddo
1219  enddo
1220 c
1221  return
1222  end
1223 c-----------------------------------------------------------------------
1224  subroutine mxmur2(a,n1,b,n2,c,n3)
1225 C----------------------------------------------------------------------
1226 C
1227 C Matrix-vector product routine.
1228 C NOTE: Use assembly coded routine if available.
1229 C
1230 C---------------------------------------------------------------------
1231  REAL A(N1,N2),B(N2,N3),C(N1,N3)
1232 C
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
1283 c
1284  subroutine mxmur2_1(a,n1,b,n2,c,n3)
1285 c
1286  real a(n1,1),b(1,n3),c(n1,n3)
1287 c
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)
1296 c
1297  real a(n1,2),b(2,n3),c(n1,n3)
1298 c
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)
1308 c
1309  real a(n1,3),b(3,n3),c(n1,n3)
1310 c
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)
1321 c
1322  real a(n1,4),b(4,n3),c(n1,n3)
1323 c
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)
1335 c
1336  real a(n1,5),b(5,n3),c(n1,n3)
1337 c
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)
1350 c
1351  real a(n1,6),b(6,n3),c(n1,n3)
1352 c
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)
1366 c
1367  real a(n1,7),b(7,n3),c(n1,n3)
1368 c
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)
1383 c
1384  real a(n1,8),b(8,n3),c(n1,n3)
1385 c
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)
1401 c
1402  real a(n1,9),b(9,n3),c(n1,n3)
1403 c
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)
1420 c
1421  real a(n1,10),b(10,n3),c(n1,n3)
1422 c
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)
1440 c
1441  real a(n1,11),b(11,n3),c(n1,n3)
1442 c
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)
1461 c
1462  real a(n1,12),b(12,n3),c(n1,n3)
1463 c
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)
1483 c
1484  real a(n1,13),b(13,n3),c(n1,n3)
1485 c
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)
1506 c
1507  real a(n1,14),b(14,n3),c(n1,n3)
1508 c
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)
1530 c
1531  real a(n1,15),b(15,n3),c(n1,n3)
1532 c
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)
1555 c
1556  real a(n1,16),b(16,n3),c(n1,n3)
1557 c
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
1580 c-----------------------------------------------------------------------
1581  subroutine mxmur3(a,n1,b,n2,c,n3)
1582 C----------------------------------------------------------------------
1583 C
1584 C Matrix-vector product routine.
1585 C NOTE: Use assembly coded routine if available.
1586 C
1587 C---------------------------------------------------------------------
1588  REAL A(N1,N2),B(N2,N3),C(N1,N3)
1589 C
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
1640 c
1641  subroutine mxmur3_16(a,n1,b,n2,c,n3)
1642  real a(n1,n2),b(n2,16),c(n1,16)
1643 c
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
1679 c
1680  enddo
1681 c
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)
1686 c
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
1720 c
1721  enddo
1722 c
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)
1727 c
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
1759 c
1760  enddo
1761 c
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)
1766 c
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
1796 c
1797  enddo
1798 c
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)
1803 c
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
1831 c
1832  enddo
1833 c
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)
1838 c
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
1864 c
1865  enddo
1866 c
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)
1871 c
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
1895 c
1896  enddo
1897 c
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)
1902 c
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
1924 c
1925  enddo
1926 c
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)
1931 c
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
1951 c
1952  enddo
1953 c
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)
1958 c
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
1976 c
1977  enddo
1978 c
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)
1983 c
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
1999 c
2000  enddo
2001 c
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)
2006 c
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
2020 c
2021  enddo
2022 c
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)
2027 c
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
2039 c
2040  enddo
2041 c
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)
2046 c
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
2056 c
2057  enddo
2058 c
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)
2063 c
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
2071 c
2072  enddo
2073 c
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)
2078 c
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
2085 c
2086  return
2087  end
2088 C----------------------------------------------------------------------
2089  subroutine mxmd(a,n1,b,n2,c,n3)
2090 C
2091 C Matrix-vector product routine.
2092 C NOTE: Use assembly coded routine if available.
2093 C
2094 C---------------------------------------------------------------------
2095  REAL A(N1,N2),B(N2,N3),C(N1,N3)
2096  REAL ONE,ZERO,EPS
2097 C
2098 C
2099 C
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
2105 c-----------------------------------------------------------------------
2106  subroutine mxmfb(a,n1,b,n2,c,n3)
2107 C-----------------------------------------------------------------------
2108 C
2109 C Matrix-vector product routine.
2110 C NOTE: Use assembly coded routine if available.
2111 C
2112 C----------------------------------------------------------------------
2113  REAL A(N1,N2),B(N2,N3),C(N1,N3)
2114 C
2115  integer wdsize
2116  save wdsize
2117  data wdsize/0/
2118 c
2119 c First call: determine word size for dgemm/sgemm discrimination, below.
2120 c
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
2127 c
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
2183 c
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
2191 c-----------------------------------------------------------------------
2192  subroutine mxmfb_1(a,n1,b,n2,c,n3)
2193 c
2194  real a(n1,1),b(1,n3),c(n1,n3)
2195 c
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
2203 c-----------------------------------------------------------------------
2204  subroutine mxmfb_2(a,n1,b,n2,c,n3)
2205 c
2206  real a(n1,2),b(2,n3),c(n1,n3)
2207 c
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
2216 c-----------------------------------------------------------------------
2217  subroutine mxmfb_3(a,n1,b,n2,c,n3)
2218 c
2219  real a(n1,3),b(3,n3),c(n1,n3)
2220 c
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
2230 c-----------------------------------------------------------------------
2231  subroutine mxmfb_4(a,n1,b,n2,c,n3)
2232 c
2233  real a(n1,4),b(4,n3),c(n1,n3)
2234 c
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
2245 c-----------------------------------------------------------------------
2246  subroutine mxmfb_5(a,n1,b,n2,c,n3)
2247 c
2248  real a(n1,5),b(5,n3),c(n1,n3)
2249 c
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
2261 c-----------------------------------------------------------------------
2262  subroutine mxmfb_6(a,n1,b,n2,c,n3)
2263 c
2264  real a(n1,6),b(6,n3),c(n1,n3)
2265 c
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
2278 c-----------------------------------------------------------------------
2279  subroutine mxmfb_7(a,n1,b,n2,c,n3)
2280 c
2281  real a(n1,7),b(7,n3),c(n1,n3)
2282 c
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
2296 c-----------------------------------------------------------------------
2297  subroutine mxmfb_8(a,n1,b,n2,c,n3)
2298 c
2299  real a(n1,8),b(8,n3),c(n1,n3)
2300 c
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
2315 c-----------------------------------------------------------------------
2316  subroutine mxmfb_9(a,n1,b,n2,c,n3)
2317 c
2318  real a(n1,9),b(9,n3),c(n1,n3)
2319 c
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
2335 c-----------------------------------------------------------------------
2336  subroutine mxmfb_10(a,n1,b,n2,c,n3)
2337 c
2338  real a(n1,10),b(10,n3),c(n1,n3)
2339 c
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
2356 c-----------------------------------------------------------------------
2357  subroutine mxmfb_11(a,n1,b,n2,c,n3)
2358 c
2359  real a(n1,11),b(11,n3),c(n1,n3)
2360 c
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
2378 c-----------------------------------------------------------------------
2379  subroutine mxmfb_12(a,n1,b,n2,c,n3)
2380 c
2381  real a(n1,12),b(12,n3),c(n1,n3)
2382 c
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
2401 c-----------------------------------------------------------------------
2402  subroutine mxmfb_13(a,n1,b,n2,c,n3)
2403 c
2404  real a(n1,13),b(13,n3),c(n1,n3)
2405 c
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
2425 c-----------------------------------------------------------------------
2426  subroutine mxmfb_14(a,n1,b,n2,c,n3)
2427 c
2428  real a(n1,14),b(14,n3),c(n1,n3)
2429 c
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
2450 c-----------------------------------------------------------------------
2451  subroutine mxmfb_15(a,n1,b,n2,c,n3)
2452 c
2453  real a(n1,15),b(15,n3),c(n1,n3)
2454 c
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
2476 c-----------------------------------------------------------------------
2477  subroutine mxmfb_16(a,n1,b,n2,c,n3)
2478 c
2479  real a(n1,16),b(16,n3),c(n1,n3)
2480 c
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
2503 c-----------------------------------------------------------------------
2504  subroutine mxmfb_17(a,n1,b,n2,c,n3)
2505 c
2506  real a(n1,17),b(17,n3),c(n1,n3)
2507 c
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
2531 c-----------------------------------------------------------------------
2532  subroutine mxmfb_18(a,n1,b,n2,c,n3)
2533 c
2534  real a(n1,18),b(18,n3),c(n1,n3)
2535 c
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
2560 c-----------------------------------------------------------------------
2561  subroutine mxmfb_19(a,n1,b,n2,c,n3)
2562 c
2563  real a(n1,19),b(19,n3),c(n1,n3)
2564 c
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
2590 c-----------------------------------------------------------------------
2591  subroutine mxmfb_20(a,n1,b,n2,c,n3)
2592 c
2593  real a(n1,20),b(20,n3),c(n1,n3)
2594 c
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
2621 c-----------------------------------------------------------------------
2622  subroutine mxmfb_21(a,n1,b,n2,c,n3)
2623 c
2624  real a(n1,21),b(21,n3),c(n1,n3)
2625 c
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
2653 c-----------------------------------------------------------------------
2654  subroutine mxmfb_22(a,n1,b,n2,c,n3)
2655 c
2656  real a(n1,22),b(22,n3),c(n1,n3)
2657 c
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
2686 c-----------------------------------------------------------------------
2687  subroutine mxmfb_23(a,n1,b,n2,c,n3)
2688 c
2689  real a(n1,23),b(23,n3),c(n1,n3)
2690 c
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
2720 c-----------------------------------------------------------------------
2721  subroutine mxmfb_24(a,n1,b,n2,c,n3)
2722 c
2723  real a(n1,24),b(24,n3),c(n1,n3)
2724 c
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
2755 c-----------------------------------------------------------------------
2756  subroutine mxmf3(a,n1,b,n2,c,n3)
2757 C-----------------------------------------------------------------------
2758 C
2759 C Matrix-vector product routine.
2760 C NOTE: Use assembly coded routine if available.
2761 C
2762 C----------------------------------------------------------------------
2763  REAL A(N1,N2),B(N2,N3),C(N1,N3)
2764 C
2765  integer wdsize
2766  save wdsize
2767  data wdsize/0/
2768 c
2769 c First call: determine word size for dgemm/sgemm discrimination, below.
2770 c
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
2777 c
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
2833 c
2834  one=1.0
2835  zero=0.0
2836  call dgemm( 'N','N',n1,n3,n2,one,a,n1,b,n2,zero,c,n1)
2837 c
2838 c N0=N1*N3
2839 c DO 10 I=1,N0
2840 c C(I,1)=0.
2841 c 10 CONTINUE
2842 c DO 100 J=1,N3
2843 c DO 100 K=1,N2
2844 c BB=B(K,J)
2845 c DO 100 I=1,N1
2846 c C(I,J)=C(I,J)+A(I,K)*BB
2847 c 100 CONTINUE
2848 
2849  endif
2850  return
2851  end
2852 c-----------------------------------------------------------------------
2853  subroutine mxmf3_1(a,n1,b,n2,c,n3)
2854 c
2855  real a(n1,1),b(1,n3),c(n1,n3)
2856 c
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
2864 c-----------------------------------------------------------------------
2865  subroutine mxmf3_2(a,n1,b,n2,c,n3)
2866 c
2867  real a(n1,2),b(2,n3),c(n1,n3)
2868 c
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
2877 c-----------------------------------------------------------------------
2878  subroutine mxmf3_3(a,n1,b,n2,c,n3)
2879 c
2880  real a(n1,3),b(3,n3),c(n1,n3)
2881 c
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
2891 c-----------------------------------------------------------------------
2892  subroutine mxmf3_4(a,n1,b,n2,c,n3)
2893 c
2894  real a(n1,4),b(4,n3),c(n1,n3)
2895 c
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
2906 c-----------------------------------------------------------------------
2907  subroutine mxmf3_5(a,n1,b,n2,c,n3)
2908 c
2909  real a(n1,5),b(5,n3),c(n1,n3)
2910 c
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
2922 c-----------------------------------------------------------------------
2923  subroutine mxmf3_6(a,n1,b,n2,c,n3)
2924 c
2925  real a(n1,6),b(6,n3),c(n1,n3)
2926 c
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
2939 c-----------------------------------------------------------------------
2940  subroutine mxmf3_7(a,n1,b,n2,c,n3)
2941 c
2942  real a(n1,7),b(7,n3),c(n1,n3)
2943 c
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
2957 c-----------------------------------------------------------------------
2958  subroutine mxmf3_8(a,n1,b,n2,c,n3)
2959 c
2960  real a(n1,8),b(8,n3),c(n1,n3)
2961 c
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
2976 c-----------------------------------------------------------------------
2977  subroutine mxmf3_9(a,n1,b,n2,c,n3)
2978 c
2979  real a(n1,9),b(9,n3),c(n1,n3)
2980 c
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
2996 c-----------------------------------------------------------------------
2997  subroutine mxmf3_10(a,n1,b,n2,c,n3)
2998 c
2999  real a(n1,10),b(10,n3),c(n1,n3)
3000 c
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
3017 c-----------------------------------------------------------------------
3018  subroutine mxmf3_11(a,n1,b,n2,c,n3)
3019 c
3020  real a(n1,11),b(11,n3),c(n1,n3)
3021 c
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
3039 c-----------------------------------------------------------------------
3040  subroutine mxmf3_12(a,n1,b,n2,c,n3)
3041 c
3042  real a(n1,12),b(12,n3),c(n1,n3)
3043 c
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
3062 c-----------------------------------------------------------------------
3063  subroutine mxmf3_13(a,n1,b,n2,c,n3)
3064 c
3065  real a(n1,13),b(13,n3),c(n1,n3)
3066 c
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
3086 c-----------------------------------------------------------------------
3087  subroutine mxmf3_14(a,n1,b,n2,c,n3)
3088 c
3089  real a(n1,14),b(14,n3),c(n1,n3)
3090 c
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
3111 c-----------------------------------------------------------------------
3112  subroutine mxmf3_15(a,n1,b,n2,c,n3)
3113 c
3114  real a(n1,15),b(15,n3),c(n1,n3)
3115 c
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
3137 c-----------------------------------------------------------------------
3138  subroutine mxmf3_16(a,n1,b,n2,c,n3)
3139 c
3140  real a(n1,16),b(16,n3),c(n1,n3)
3141 c
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
3164 c-----------------------------------------------------------------------
3165  subroutine mxmf3_17(a,n1,b,n2,c,n3)
3166 c
3167  real a(n1,17),b(17,n3),c(n1,n3)
3168 c
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
3192 c-----------------------------------------------------------------------
3193  subroutine mxmf3_18(a,n1,b,n2,c,n3)
3194 c
3195  real a(n1,18),b(18,n3),c(n1,n3)
3196 c
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
3221 c-----------------------------------------------------------------------
3222  subroutine mxmf3_19(a,n1,b,n2,c,n3)
3223 c
3224  real a(n1,19),b(19,n3),c(n1,n3)
3225 c
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
3251 c-----------------------------------------------------------------------
3252  subroutine mxmf3_20(a,n1,b,n2,c,n3)
3253 c
3254  real a(n1,20),b(20,n3),c(n1,n3)
3255 c
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
3282 c-----------------------------------------------------------------------
3283  subroutine mxmf3_21(a,n1,b,n2,c,n3)
3284 c
3285  real a(n1,21),b(21,n3),c(n1,n3)
3286 c
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
3314 c-----------------------------------------------------------------------
3315  subroutine mxmf3_22(a,n1,b,n2,c,n3)
3316 c
3317  real a(n1,22),b(22,n3),c(n1,n3)
3318 c
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
3347 c-----------------------------------------------------------------------
3348  subroutine mxmf3_23(a,n1,b,n2,c,n3)
3349 c
3350  real a(n1,23),b(23,n3),c(n1,n3)
3351 c
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
3381 c-----------------------------------------------------------------------
3382  subroutine mxmf3_24(a,n1,b,n2,c,n3)
3383 c
3384  real a(n1,24),b(24,n3),c(n1,n3)
3385 c
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
3416 c-----------------------------------------------------------------------
3417  subroutine mxm44(a,n1,b,n2,c,n3)
3418 C-----------------------------------------------------------------------
3419 C
3420 C NOTE -- this code has been set up with the "mxmf3" routine
3421 c referenced in memtime.f. On most machines, the f2
3422 c and f3 versions give the same performance (f2 is the
3423 c nekton standard). On the t3e, f3 is noticeably faster.
3424 c pff 10/5/98
3425 C
3426 C
3427 C Matrix-vector product routine.
3428 C NOTE: Use assembly coded routine if available.
3429 C
3430 C----------------------------------------------------------------------
3431  REAL A(N1,N2),B(N2,N3),C(N1,N3)
3432 c
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
3440 c
3441  return
3442  end
3443 c
3444 c-----------------------------------------------------------------------
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
3451 c
3452 c matrix multiply with a 4x4 pencil
3453 c
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
3875 c-----------------------------------------------------------------------
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
3921 c-----------------------------------------------------------------------
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 
3950 c 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 
3971 c-------------------------------------------------------
3972 c mem test
3973 c-------------------------------------------------------
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
3987 c write(6,*) loop,time,time2,time1,overh
3988  flops=iops/(1.0e6*time)
3989  s(n,1) = flops
3990 c
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')
3994 c
3995 c-------------------------------------------------------
3996 c fast test
3997 c-------------------------------------------------------
3998 c
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
4013 c
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')
4016 c
4017  endif
4018  enddo
4019 c
4020  return
4021  end
4022 c-----------------------------------------------------------------------
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)
4029 c ^ ^ ^ |__ N^2xN, NxN, NxN^2
4030 c matrix order N __| | |__________which mxm
4031 c |
4032 c |__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)
4067 c
4068 c Case by case
4069 c
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
4088 c-----------------------------------------------------------------------
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
4116 c-----------------------------------------------------------------------
4117  subroutine mxma(a,n1,b,n2,c,n3)
4118 C
4119 C Compute C = A*B , for contiguously packed matrices A,B, and C.
4120 C
4121 C
4122 C Compile with -r8 option for 64-bit arithmetic, or convert real
4123 c to real*8
4124 c
4125  real a(n1,n2),b(n2,n3),c(n1,n3)
4126 c
4127  include 'SIZE'
4128  include 'OPCTR'
4129  include 'TOTAL'
4130 c
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
4143 c
4144 c
4145  call mxma2(a,n1,b,n2,c,n3) ! In some cases, this is faster.
4146 c
4147  return
4148  end
4149 c-----------------------------------------------------------------------
4150  subroutine mxma2(a,n1,b,n2,c,n3)
4151 c
4152 c unrolled loop routine written by paul fischer
4153 c
4154 c this version computines C = C + A x B
4155 c
4156  real a(n1,n2),b(n2,n3),c(n1,n3)
4157 C
4158  integer wdsize
4159  save wdsize
4160  data wdsize/0/
4161 c
4162 c First call: determine word size for dgemm/sgemm discrimination, below.
4163 c
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
4170 c
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
4228 c
4229  return
4230  end
4231 c-----------------------------------------------------------------------
4232  subroutine mxa1(a,n1,b,n2,c,n3)
4233 c
4234  real a(n1,1),b(1,n3),c(n1,n3)
4235 c
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
4244 c-----------------------------------------------------------------------
4245  subroutine mxa2(a,n1,b,n2,c,n3)
4246 c
4247  real a(n1,2),b(2,n3),c(n1,n3)
4248 c
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
4258 c-----------------------------------------------------------------------
4259  subroutine mxa3(a,n1,b,n2,c,n3)
4260 c
4261  real a(n1,3),b(3,n3),c(n1,n3)
4262 c
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
4273 c-----------------------------------------------------------------------
4274  subroutine mxa4(a,n1,b,n2,c,n3)
4275 c
4276  real a(n1,4),b(4,n3),c(n1,n3)
4277 c
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
4289 c-----------------------------------------------------------------------
4290  subroutine mxa5(a,n1,b,n2,c,n3)
4291 c
4292  real a(n1,5),b(5,n3),c(n1,n3)
4293 c
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
4306 c-----------------------------------------------------------------------
4307  subroutine mxa6(a,n1,b,n2,c,n3)
4308 c
4309  real a(n1,6),b(6,n3),c(n1,n3)
4310 c
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
4324 c-----------------------------------------------------------------------
4325  subroutine mxa7(a,n1,b,n2,c,n3)
4326 c
4327  real a(n1,7),b(7,n3),c(n1,n3)
4328 c
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
4343 c-----------------------------------------------------------------------
4344  subroutine mxa8(a,n1,b,n2,c,n3)
4345 c
4346  real a(n1,8),b(8,n3),c(n1,n3)
4347 c
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
4363 c-----------------------------------------------------------------------
4364  subroutine mxa9(a,n1,b,n2,c,n3)
4365 c
4366  real a(n1,9),b(9,n3),c(n1,n3)
4367 c
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
4384 c-----------------------------------------------------------------------
4385  subroutine mxa10(a,n1,b,n2,c,n3)
4386 c
4387  real a(n1,10),b(10,n3),c(n1,n3)
4388 c
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
4406 c-----------------------------------------------------------------------
4407  subroutine mxa11(a,n1,b,n2,c,n3)
4408 c
4409  real a(n1,11),b(11,n3),c(n1,n3)
4410 c
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
4429 c-----------------------------------------------------------------------
4430  subroutine mxa12(a,n1,b,n2,c,n3)
4431 c
4432  real a(n1,12),b(12,n3),c(n1,n3)
4433 c
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
4453 c-----------------------------------------------------------------------
4454  subroutine mxa13(a,n1,b,n2,c,n3)
4455 c
4456  real a(n1,13),b(13,n3),c(n1,n3)
4457 c
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
4478 c-----------------------------------------------------------------------
4479  subroutine mxa14(a,n1,b,n2,c,n3)
4480 c
4481  real a(n1,14),b(14,n3),c(n1,n3)
4482 c
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
4504 c-----------------------------------------------------------------------
4505  subroutine mxa15(a,n1,b,n2,c,n3)
4506 c
4507  real a(n1,15),b(15,n3),c(n1,n3)
4508 c
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
4531 c-----------------------------------------------------------------------
4532  subroutine mxa16(a,n1,b,n2,c,n3)
4533 c
4534  real a(n1,16),b(16,n3),c(n1,n3)
4535 c
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
4559 c-----------------------------------------------------------------------
4560  subroutine mxa17(a,n1,b,n2,c,n3)
4561 c
4562  real a(n1,17),b(17,n3),c(n1,n3)
4563 c
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
4588 c-----------------------------------------------------------------------
4589  subroutine mxa18(a,n1,b,n2,c,n3)
4590 c
4591  real a(n1,18),b(18,n3),c(n1,n3)
4592 c
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
4618 c-----------------------------------------------------------------------
4619  subroutine mxa19(a,n1,b,n2,c,n3)
4620 c
4621  real a(n1,19),b(19,n3),c(n1,n3)
4622 c
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
4649 c-----------------------------------------------------------------------
4650  subroutine mxa20(a,n1,b,n2,c,n3)
4651 c
4652  real a(n1,20),b(20,n3),c(n1,n3)
4653 c
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
4681 c-----------------------------------------------------------------------
4682  subroutine mxa21(a,n1,b,n2,c,n3)
4683 c
4684  real a(n1,21),b(21,n3),c(n1,n3)
4685 c
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
4714 c-----------------------------------------------------------------------
4715  subroutine mxa22(a,n1,b,n2,c,n3)
4716 c
4717  real a(n1,22),b(22,n3),c(n1,n3)
4718 c
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
4748 c-----------------------------------------------------------------------
4749  subroutine mxa23(a,n1,b,n2,c,n3)
4750 c
4751  real a(n1,23),b(23,n3),c(n1,n3)
4752 c
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
4783 c-----------------------------------------------------------------------
4784  subroutine mxa24(a,n1,b,n2,c,n3)
4785 c
4786  real a(n1,24),b(24,n3),c(n1,n3)
4787 c
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
4819 c-----------------------------------------------------------------------
4820  subroutine mxm44_0a(a, m, b, k, c, n)
4821 c
4822 c routine written by Bruce Curtiss
4823 c
4824 c subroutine matmul44(m, n, k, a, lda, b, ldb, c, ldc)
4825 c real a(lda,k), b(ldb,n), c(ldc,n)
4826 c
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
4830 c
4831 c matrix multiply with a 4x4 pencil
4832 c
4833 
4834  mresid = iand(m,3)
4835  nresid = iand(n,3)
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)
4861 c
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
4904 c 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 
4997 c 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
5017 c 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
5081 c 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
5134 c 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
5182 c 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
5256 c-----------------------------------------------------------------------
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
5312 c-----------------------------------------------------------------------
nekgsync
subroutine nekgsync()
Definition: comm_mpi.f:502
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dgemm.f:3
out_anal
subroutine out_anal(s, a, nn, c, nt, itmax, name8, k, ivb)
Definition: mxm_std.f:4090
mxf15
subroutine mxf15(a, n1, b, n2, c, n3)
Definition: mxm_std.f:329
mxm44_2_t
subroutine mxm44_2_t(a, m, b, k, c, n)
Definition: mxm_std.f:3877
mxmur2_6
subroutine mxmur2_6(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1350
mxmfb_17
subroutine mxmfb_17(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2505
mxa24
subroutine mxa24(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4785
mxmur2
subroutine mxmur2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1225
mxmfb_9
subroutine mxmfb_9(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2317
madd
subroutine madd(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1212
mxmf3_6
subroutine mxmf3_6(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2924
mxmfb_13
subroutine mxmfb_13(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2403
mxa22
subroutine mxa22(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4716
mxf7
subroutine mxf7(a, n1, b, n2, c, n3)
Definition: mxm_std.f:157
mxmf3_11
subroutine mxmf3_11(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3019
mxmur2_10
subroutine mxmur2_10(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1420
mxa7
subroutine mxa7(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4326
mxa5
subroutine mxa5(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4291
mxf22
subroutine mxf22(a, n1, b, n2, c, n3)
Definition: mxm_std.f:532
mxmf3_5
subroutine mxmf3_5(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2908
mxmfb_19
subroutine mxmfb_19(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2562
mxa14
subroutine mxa14(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4480
mxf11
subroutine mxf11(a, n1, b, n2, c, n3)
Definition: mxm_std.f:235
mxm44_0a
subroutine mxm44_0a(a, m, b, k, c, n)
Definition: mxm_std.f:4821
mxa12
subroutine mxa12(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4431
mxa11
subroutine mxa11(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4408
mxmf3_21
subroutine mxmf3_21(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3284
mxa17
subroutine mxa17(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4561
mxmur2_15
subroutine mxmur2_15(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1530
mxf19
subroutine mxf19(a, n1, b, n2, c, n3)
Definition: mxm_std.f:439
mxmur3
subroutine mxmur3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1582
mxmur3_7
subroutine mxmur3_7(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1957
mxmfb_15
subroutine mxmfb_15(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2452
mxmf3_14
subroutine mxmf3_14(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3088
mxmur3_11
subroutine mxmur3_11(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1837
mxmf3_19
subroutine mxmf3_19(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3223
mxmur2_16
subroutine mxmur2_16(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1555
mxa9
subroutine mxa9(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4365
mxmf3_9
subroutine mxmf3_9(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2978
mxmfb_8
subroutine mxmfb_8(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2298
mxm_test_all
subroutine mxm_test_all(nid, ivb)
Definition: mxm_std.f:1109
mxf4
subroutine mxf4(a, n1, b, n2, c, n3)
Definition: mxm_std.f:109
mxf16
subroutine mxf16(a, n1, b, n2, c, n3)
Definition: mxm_std.f:355
mxmf3_1
subroutine mxmf3_1(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2854
mxmur3_16
subroutine mxmur3_16(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1642
mxf17
subroutine mxf17(a, n1, b, n2, c, n3)
Definition: mxm_std.f:382
mxf6
subroutine mxf6(a, n1, b, n2, c, n3)
Definition: mxm_std.f:140
mxmfb_11
subroutine mxmfb_11(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2358
mxmf3_3
subroutine mxmf3_3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2879
mxf23
subroutine mxf23(a, n1, b, n2, c, n3)
Definition: mxm_std.f:565
mxa23
subroutine mxa23(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4750
mxmur2_12
subroutine mxmur2_12(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1461
mxmur2_5
subroutine mxmur2_5(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1335
mxmf3_18
subroutine mxmf3_18(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3194
mxmfb_5
subroutine mxmfb_5(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2247
mxmur3_9
subroutine mxmur3_9(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1901
mxf18
subroutine mxf18(a, n1, b, n2, c, n3)
Definition: mxm_std.f:410
mxa6
subroutine mxa6(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4308
mxmfb_1
subroutine mxmfb_1(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2193
mxmur2_8
subroutine mxmur2_8(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1383
mxmfb
subroutine mxmfb(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2107
mxmf3_22
subroutine mxmf3_22(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3316
mxf10
subroutine mxf10(a, n1, b, n2, c, n3)
Definition: mxm_std.f:214
mxmf3_15
subroutine mxmf3_15(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3113
mxf5
subroutine mxf5(a, n1, b, n2, c, n3)
Definition: mxm_std.f:124
mxa3
subroutine mxa3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4260
mxa20
subroutine mxa20(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4651
mxmfb_16
subroutine mxmfb_16(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2478
mxa2
subroutine mxa2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4246
mxa8
subroutine mxa8(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4345
mxm44_0_t
subroutine mxm44_0_t(a, m, b, k, c, n)
Definition: mxm_std.f:3446
mxf13
subroutine mxf13(a, n1, b, n2, c, n3)
Definition: mxm_std.f:280
mxf14
subroutine mxf14(a, n1, b, n2, c, n3)
Definition: mxm_std.f:304
mxmur3_15
subroutine mxmur3_15(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1685
mxa16
subroutine mxa16(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4533
mxmf3_8
subroutine mxmf3_8(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2959
mxmfb_22
subroutine mxmfb_22(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2655
mxmur3_10
subroutine mxmur3_10(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1870
mxmur3_6
subroutine mxmur3_6(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1982
mxmfb_6
subroutine mxmfb_6(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2263
mxmf3_23
subroutine mxmf3_23(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3349
mxf1
subroutine mxf1(a, n1, b, n2, c, n3)
Definition: mxm_std.f:70
mxa1
subroutine mxa1(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4233
mxmur2_2
subroutine mxmur2_2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1296
mxf21
subroutine mxf21(a, n1, b, n2, c, n3)
Definition: mxm_std.f:500
mxa10
subroutine mxa10(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4386
mxmur2_1
subroutine mxmur2_1(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1285
mxmfb_3
subroutine mxmfb_3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2218
mxm44
subroutine mxm44(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3418
mxm_analyze
subroutine mxm_analyze(s, a, nn, c, nt, ivb)
Definition: mxm_std.f:4024
mxmur3_5
subroutine mxmur3_5(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2005
mxa18
subroutine mxa18(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4590
mxf9
subroutine mxf9(a, n1, b, n2, c, n3)
Definition: mxm_std.f:194
mxmfb_23
subroutine mxmfb_23(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2688
mxmf3_7
subroutine mxmf3_7(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2941
mxmf3_10
subroutine mxmf3_10(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2998
mxmur2_14
subroutine mxmur2_14(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1506
mxmur2_13
subroutine mxmur2_13(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1483
mxmf3_4
subroutine mxmf3_4(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2893
mxmu4
subroutine mxmu4(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1182
mxmur2_4
subroutine mxmur2_4(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1321
mxmfb_24
subroutine mxmfb_24(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2722
mxmf3_12
subroutine mxmf3_12(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3041
mxmfb_18
subroutine mxmfb_18(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2533
mxmf2
subroutine mxmf2(A, N1, B, N2, C, N3)
Definition: mxm_std.f:3
mxmtest
subroutine mxmtest(s, nn, cn, mxmt, name, k, ivb)
Definition: mxm_std.f:3923
mxa21
subroutine mxa21(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4683
mxmfb_10
subroutine mxmfb_10(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2337
mxms
subroutine mxms(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1160
mxmf3
subroutine mxmf3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2757
mxmfb_12
subroutine mxmfb_12(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2380
mxm44_2a
subroutine mxm44_2a(a, m, b, k, c, n)
Definition: mxm_std.f:5258
mxf3
subroutine mxf3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:95
mxmur3_14
subroutine mxmur3_14(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1726
initab
subroutine initab(a, b, n)
Definition: mxm_std.f:1145
mxmf3_17
subroutine mxmf3_17(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3166
mxm44_0
subroutine mxm44_0(a, m, b, k, c, n)
Definition: mxm_std.f:634
mxa19
subroutine mxa19(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4620
mxf2
subroutine mxf2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:82
mxmur3_2
subroutine mxmur3_2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2062
mxmur2_9
subroutine mxmur2_9(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1401
mxmfb_20
subroutine mxmfb_20(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2592
mxmur3_13
subroutine mxmur3_13(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1765
mxmf3_20
subroutine mxmf3_20(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3253
mxmf3_13
subroutine mxmf3_13(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3064
mxmur3_1
subroutine mxmur3_1(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2077
mxa15
subroutine mxa15(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4506
mxmfb_2
subroutine mxmfb_2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2205
mxmur3_12
subroutine mxmur3_12(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1802
mxmfb_21
subroutine mxmfb_21(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2623
mxmur2_3
subroutine mxmur2_3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1308
mxma
subroutine mxma(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4118
mxmur3_8
subroutine mxmur3_8(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1930
mxmf3_2
subroutine mxmf3_2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2866
mxmfb_14
subroutine mxmfb_14(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2427
mxmfb_7
subroutine mxmfb_7(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2280
mxmf3_24
subroutine mxmf3_24(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3383
mxf24
subroutine mxf24(a, n1, b, n2, c, n3)
Definition: mxm_std.f:599
mxm44_2
subroutine mxm44_2(a, m, b, k, c, n)
Definition: mxm_std.f:1063
mxma2
subroutine mxma2(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4151
mxmur2_7
subroutine mxmur2_7(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1366
mxmd
subroutine mxmd(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2090
mxf8
subroutine mxf8(a, n1, b, n2, c, n3)
Definition: mxm_std.f:175
mxmur3_3
subroutine mxmur3_3(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2045
mxa4
subroutine mxa4(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4275
mxf20
subroutine mxf20(a, n1, b, n2, c, n3)
Definition: mxm_std.f:469
mxa13
subroutine mxa13(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4455
mxmf3_16
subroutine mxmf3_16(a, n1, b, n2, c, n3)
Definition: mxm_std.f:3139
mxf12
subroutine mxf12(a, n1, b, n2, c, n3)
Definition: mxm_std.f:257
mxmfb_4
subroutine mxmfb_4(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2232
mxmur3_4
subroutine mxmur3_4(a, n1, b, n2, c, n3)
Definition: mxm_std.f:2026
mxmur2_11
subroutine mxmur2_11(a, n1, b, n2, c, n3)
Definition: mxm_std.f:1440
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2