Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
mxm_std.F90
Go to the documentation of this file.
1 
2 subroutine mxmf2(A,N1,B,N2,C,N3)
3  use kinds, only : dp
4  integer :: n1, n2, n3
5  real(DP) :: a(n1,n2),b(n2,n3),c(n1,n3)
6 
7  select case (n2)
8  case (1 : 8)
9  select case (n2)
10  case (8)
11  call mxf8(a,n1,b,n2,c,n3)
12  case (1)
13  call mxf1(a,n1,b,n2,c,n3)
14  case (2)
15  call mxf2(a,n1,b,n2,c,n3)
16  case (3)
17  call mxf3(a,n1,b,n2,c,n3)
18  case (4)
19  call mxf4(a,n1,b,n2,c,n3)
20  case (5)
21  call mxf5(a,n1,b,n2,c,n3)
22  case (6)
23  call mxf6(a,n1,b,n2,c,n3)
24  case (7)
25  call mxf7(a,n1,b,n2,c,n3)
26  end select
27  case (9 : 16)
28  select case (n2)
29  case (12)
30  call mxf12(a,n1,b,n2,c,n3)
31  case (9)
32  call mxf9(a,n1,b,n2,c,n3)
33  case (10)
34  call mxf10(a,n1,b,n2,c,n3)
35  case (11)
36  call mxf11(a,n1,b,n2,c,n3)
37  case (13)
38  call mxf13(a,n1,b,n2,c,n3)
39  case (14)
40  call mxf14(a,n1,b,n2,c,n3)
41  case (15)
42  call mxf15(a,n1,b,n2,c,n3)
43  case (16)
44  call mxf16(a,n1,b,n2,c,n3)
45  end select
46  case (17 : 24)
47  select case (n2)
48  case (17)
49  call mxf17(a,n1,b,n2,c,n3)
50  case (18)
51  call mxf18(a,n1,b,n2,c,n3)
52  case (19)
53  call mxf19(a,n1,b,n2,c,n3)
54  case (20)
55  call mxf20(a,n1,b,n2,c,n3)
56  case (21)
57  call mxf21(a,n1,b,n2,c,n3)
58  case (22)
59  call mxf22(a,n1,b,n2,c,n3)
60  case (23)
61  call mxf23(a,n1,b,n2,c,n3)
62  case (24)
63  call mxf24(a,n1,b,n2,c,n3)
64  end select
65  case default
66  call mxm44_0(a,n1,b,n2,c,n3)
67  end select
68 
69  return
70 end subroutine mxmf2
71 
72 !-----------------------------------------------------------------------
73  subroutine mxf1(a,n1,b,n2,c,n3)
74 
75  use kinds, only : dp
76  real(DP) :: a(n1,1),b(1,n3),c(n1,n3)
77 
78  do j=1,n3
79  do i=1,n1
80  c(i,j) = a(i,1)*b(1,j)
81  enddo
82  enddo
83  return
84  end subroutine mxf1
85 !-----------------------------------------------------------------------
86  subroutine mxf2(a,n1,b,n2,c,n3)
87 
88  use kinds, only : dp
89  real(DP) :: a(n1,2),b(2,n3),c(n1,n3)
90 
91  do j=1,n3
92  do i=1,n1
93  c(i,j) = a(i,1)*b(1,j) &
94  + a(i,2)*b(2,j)
95  enddo
96  enddo
97  return
98  end subroutine mxf2
99 !-----------------------------------------------------------------------
100  subroutine mxf3(a,n1,b,n2,c,n3)
101 
102  use kinds, only : dp
103  real(DP) :: a(n1,3),b(3,n3),c(n1,n3)
104 
105  do j=1,n3
106  do i=1,n1
107  c(i,j) = a(i,1)*b(1,j) &
108  + a(i,2)*b(2,j) &
109  + a(i,3)*b(3,j)
110  enddo
111  enddo
112  return
113  end subroutine mxf3
114 !-----------------------------------------------------------------------
115  subroutine mxf4(a,n1,b,n2,c,n3)
116 
117  use kinds, only : dp
118  real(DP) :: a(n1,4),b(4,n3),c(n1,n3)
119 
120  do j=1,n3
121  do i=1,n1
122  c(i,j) = a(i,1)*b(1,j) &
123  + a(i,2)*b(2,j) &
124  + a(i,3)*b(3,j) &
125  + a(i,4)*b(4,j)
126  enddo
127  enddo
128  return
129  end subroutine mxf4
130 !-----------------------------------------------------------------------
131  subroutine mxf5(a,n1,b,n2,c,n3)
132 
133  use kinds, only : dp
134  real(DP) :: a(n1,5),b(5,n3),c(n1,n3)
135 
136  do j=1,n3
137  do i=1,n1
138  c(i,j) = a(i,1)*b(1,j) &
139  + a(i,2)*b(2,j) &
140  + a(i,3)*b(3,j) &
141  + a(i,4)*b(4,j) &
142  + a(i,5)*b(5,j)
143  enddo
144  enddo
145  return
146  end subroutine mxf5
147 !-----------------------------------------------------------------------
148  subroutine mxf6(a,n1,b,n2,c,n3)
149 
150  use kinds, only : dp
151  real(DP) :: a(n1,6),b(6,n3),c(n1,n3)
152 
153  do j=1,n3
154  do i=1,n1
155  c(i,j) = a(i,1)*b(1,j) &
156  + a(i,2)*b(2,j) &
157  + a(i,3)*b(3,j) &
158  + a(i,4)*b(4,j) &
159  + a(i,5)*b(5,j) &
160  + a(i,6)*b(6,j)
161  enddo
162  enddo
163  return
164  end subroutine mxf6
165 !-----------------------------------------------------------------------
166  subroutine mxf7(a,n1,b,n2,c,n3)
167 
168  use kinds, only : dp
169  real(DP) :: a(n1,7),b(7,n3),c(n1,n3)
170 
171  do j=1,n3
172  do i=1,n1
173  c(i,j) = a(i,1)*b(1,j) &
174  + a(i,2)*b(2,j) &
175  + a(i,3)*b(3,j) &
176  + a(i,4)*b(4,j) &
177  + a(i,5)*b(5,j) &
178  + a(i,6)*b(6,j) &
179  + a(i,7)*b(7,j)
180  enddo
181  enddo
182  return
183  end subroutine mxf7
184 !-----------------------------------------------------------------------
185  subroutine mxf8(a,n1,b,n2,c,n3)
186 
187  use kinds, only : dp
188  real(DP) :: a(n1,8),b(8,n3),c(n1,n3)
189 
190  do j=1,n3
191  do i=1,n1
192  c(i,j) = a(i,1)*b(1,j) &
193  + a(i,2)*b(2,j) &
194  + a(i,3)*b(3,j) &
195  + a(i,4)*b(4,j) &
196  + a(i,5)*b(5,j) &
197  + a(i,6)*b(6,j) &
198  + a(i,7)*b(7,j) &
199  + a(i,8)*b(8,j)
200  enddo
201  enddo
202  return
203  end subroutine mxf8
204 !-----------------------------------------------------------------------
205  subroutine mxf9(a,n1,b,n2,c,n3)
206 
207  use kinds, only : dp
208  real(DP) :: a(n1,9),b(9,n3),c(n1,n3)
209 
210  do j=1,n3
211  do i=1,n1
212  c(i,j) = a(i,1)*b(1,j) &
213  + a(i,2)*b(2,j) &
214  + a(i,3)*b(3,j) &
215  + a(i,4)*b(4,j) &
216  + a(i,5)*b(5,j) &
217  + a(i,6)*b(6,j) &
218  + a(i,7)*b(7,j) &
219  + a(i,8)*b(8,j) &
220  + a(i,9)*b(9,j)
221  enddo
222  enddo
223  return
224  end subroutine mxf9
225 !-----------------------------------------------------------------------
226  subroutine mxf10(a,n1,b,n2,c,n3)
227 
228  use kinds, only : dp
229  real(DP) :: a(n1,10),b(10,n3),c(n1,n3)
230 
231  do j=1,n3
232  do i=1,n1
233  c(i,j) = a(i,1)*b(1,j) &
234  + a(i,2)*b(2,j) &
235  + a(i,3)*b(3,j) &
236  + a(i,4)*b(4,j) &
237  + a(i,5)*b(5,j) &
238  + a(i,6)*b(6,j) &
239  + a(i,7)*b(7,j) &
240  + a(i,8)*b(8,j) &
241  + a(i,9)*b(9,j) &
242  + a(i,10)*b(10,j)
243  enddo
244  enddo
245  return
246  end subroutine mxf10
247 !-----------------------------------------------------------------------
248  subroutine mxf11(a,n1,b,n2,c,n3)
249 
250  use kinds, only : dp
251  real(DP) :: a(n1,11),b(11,n3),c(n1,n3)
252 
253  do j=1,n3
254  do i=1,n1
255  c(i,j) = a(i,1)*b(1,j) &
256  + a(i,2)*b(2,j) &
257  + a(i,3)*b(3,j) &
258  + a(i,4)*b(4,j) &
259  + a(i,5)*b(5,j) &
260  + a(i,6)*b(6,j) &
261  + a(i,7)*b(7,j) &
262  + a(i,8)*b(8,j) &
263  + a(i,9)*b(9,j) &
264  + a(i,10)*b(10,j) &
265  + a(i,11)*b(11,j)
266  enddo
267  enddo
268  return
269  end subroutine mxf11
270 !-----------------------------------------------------------------------
271  subroutine mxf12(a,n1,b,n2,c,n3)
272 
273  use kinds, only : dp
274  real(DP) :: a(n1,12),b(12,n3),c(n1,n3)
275 
276  do j=1,n3
277  do i=1,n1
278  c(i,j) = a(i,1)*b(1,j) &
279  + a(i,2)*b(2,j) &
280  + a(i,3)*b(3,j) &
281  + a(i,4)*b(4,j) &
282  + a(i,5)*b(5,j) &
283  + a(i,6)*b(6,j) &
284  + a(i,7)*b(7,j) &
285  + a(i,8)*b(8,j) &
286  + a(i,9)*b(9,j) &
287  + a(i,10)*b(10,j) &
288  + a(i,11)*b(11,j) &
289  + a(i,12)*b(12,j)
290  enddo
291  enddo
292  return
293  end subroutine mxf12
294 !-----------------------------------------------------------------------
295  subroutine mxf13(a,n1,b,n2,c,n3)
296 
297  use kinds, only : dp
298  real(DP) :: a(n1,13),b(13,n3),c(n1,n3)
299 
300  do j=1,n3
301  do i=1,n1
302  c(i,j) = a(i,1)*b(1,j) &
303  + a(i,2)*b(2,j) &
304  + a(i,3)*b(3,j) &
305  + a(i,4)*b(4,j) &
306  + a(i,5)*b(5,j) &
307  + a(i,6)*b(6,j) &
308  + a(i,7)*b(7,j) &
309  + a(i,8)*b(8,j) &
310  + a(i,9)*b(9,j) &
311  + a(i,10)*b(10,j) &
312  + a(i,11)*b(11,j) &
313  + a(i,12)*b(12,j) &
314  + a(i,13)*b(13,j)
315  enddo
316  enddo
317  return
318  end subroutine mxf13
319 !-----------------------------------------------------------------------
320  subroutine mxf14(a,n1,b,n2,c,n3)
321 use kinds, only : dp
322 
323  real(DP) :: a(n1,14),b(14,n3),c(n1,n3)
324 
325  do j=1,n3
326  do i=1,n1
327  c(i,j) = a(i,1)*b(1,j) &
328  + a(i,2)*b(2,j) &
329  + a(i,3)*b(3,j) &
330  + a(i,4)*b(4,j) &
331  + a(i,5)*b(5,j) &
332  + a(i,6)*b(6,j) &
333  + a(i,7)*b(7,j) &
334  + a(i,8)*b(8,j) &
335  + a(i,9)*b(9,j) &
336  + a(i,10)*b(10,j) &
337  + a(i,11)*b(11,j) &
338  + a(i,12)*b(12,j) &
339  + a(i,13)*b(13,j) &
340  + a(i,14)*b(14,j)
341  enddo
342  enddo
343  return
344  end subroutine mxf14
345 !-----------------------------------------------------------------------
346  subroutine mxf15(a,n1,b,n2,c,n3)
347 
348  use kinds, only : dp
349  real(DP) :: a(n1,15),b(15,n3),c(n1,n3)
350 
351  do j=1,n3
352  do i=1,n1
353  c(i,j) = a(i,1)*b(1,j) &
354  + a(i,2)*b(2,j) &
355  + a(i,3)*b(3,j) &
356  + a(i,4)*b(4,j) &
357  + a(i,5)*b(5,j) &
358  + a(i,6)*b(6,j) &
359  + a(i,7)*b(7,j) &
360  + a(i,8)*b(8,j) &
361  + a(i,9)*b(9,j) &
362  + a(i,10)*b(10,j) &
363  + a(i,11)*b(11,j) &
364  + a(i,12)*b(12,j) &
365  + a(i,13)*b(13,j) &
366  + a(i,14)*b(14,j) &
367  + a(i,15)*b(15,j)
368  enddo
369  enddo
370  return
371  end subroutine mxf15
372 !-----------------------------------------------------------------------
373  subroutine mxf16(a,n1,b,n2,c,n3)
374 
375  use kinds, only : dp
376  real(DP) :: a(n1,16),b(16,n3),c(n1,n3)
377 
378  do j=1,n3
379  do i=1,n1
380 #if 1
381  c(i,j) = 0
382  do k=1,16
383  c(i,j) = c(i,j) + a(i,k)*b(k,j)
384  enddo
385 #else
386  c(i,j) = a(i,1)*b(1,j) &
387  + a(i,2)*b(2,j) &
388  + a(i,3)*b(3,j) &
389  + a(i,4)*b(4,j) &
390  + a(i,5)*b(5,j) &
391  + a(i,6)*b(6,j) &
392  + a(i,7)*b(7,j) &
393  + a(i,8)*b(8,j) &
394  + a(i,9)*b(9,j) &
395  + a(i,10)*b(10,j) &
396  + a(i,11)*b(11,j) &
397  + a(i,12)*b(12,j) &
398  + a(i,13)*b(13,j) &
399  + a(i,14)*b(14,j) &
400  + a(i,15)*b(15,j) &
401  + a(i,16)*b(16,j)
402 #endif
403  enddo
404  enddo
405  return
406  end subroutine mxf16
407 !-----------------------------------------------------------------------
408  subroutine mxf17(a,n1,b,n2,c,n3)
409 
410  use kinds, only : dp
411  real(DP) :: a(n1,17),b(17,n3),c(n1,n3)
412 
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  enddo
433  enddo
434  return
435  end subroutine mxf17
436 !-----------------------------------------------------------------------
437  subroutine mxf18(a,n1,b,n2,c,n3)
438 
439  use kinds, only : dp
440  real(DP) :: a(n1,18),b(18,n3),c(n1,n3)
441 
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  enddo
463  enddo
464  return
465  end subroutine mxf18
466 !-----------------------------------------------------------------------
467  subroutine mxf19(a,n1,b,n2,c,n3)
468 
469  use kinds, only : dp
470  real(DP) :: a(n1,19),b(19,n3),c(n1,n3)
471 
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  enddo
494  enddo
495  return
496  end subroutine mxf19
497 !-----------------------------------------------------------------------
498  subroutine mxf20(a,n1,b,n2,c,n3)
499 
500  use kinds, only : dp
501  real(DP) :: a(n1,20),b(20,n3),c(n1,n3)
502 
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  enddo
526  enddo
527  return
528  end subroutine mxf20
529 !-----------------------------------------------------------------------
530  subroutine mxf21(a,n1,b,n2,c,n3)
531 
532  use kinds, only : dp
533  real(DP) :: a(n1,21),b(21,n3),c(n1,n3)
534 
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  enddo
559  enddo
560  return
561  end subroutine mxf21
562 !-----------------------------------------------------------------------
563  subroutine mxf22(a,n1,b,n2,c,n3)
564 
565  use kinds, only : dp
566  real(DP) :: a(n1,22),b(22,n3),c(n1,n3)
567 
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  enddo
593  enddo
594  return
595  end subroutine mxf22
596 !-----------------------------------------------------------------------
597  subroutine mxf23(a,n1,b,n2,c,n3)
598 
599  use kinds, only : dp
600  real(DP) :: a(n1,23),b(23,n3),c(n1,n3)
601 
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  enddo
628  enddo
629  return
630  end subroutine mxf23
631 !-----------------------------------------------------------------------
632  subroutine mxf24(a,n1,b,n2,c,n3)
633 
634  use kinds, only : dp
635  real(DP) :: a(n1,24),b(24,n3),c(n1,n3)
636 
637  do j=1,n3
638  do i=1,n1
639  c(i,j) = a(i,1)*b(1,j) &
640  + a(i,2)*b(2,j) &
641  + a(i,3)*b(3,j) &
642  + a(i,4)*b(4,j) &
643  + a(i,5)*b(5,j) &
644  + a(i,6)*b(6,j) &
645  + a(i,7)*b(7,j) &
646  + a(i,8)*b(8,j) &
647  + a(i,9)*b(9,j) &
648  + a(i,10)*b(10,j) &
649  + a(i,11)*b(11,j) &
650  + a(i,12)*b(12,j) &
651  + a(i,13)*b(13,j) &
652  + a(i,14)*b(14,j) &
653  + a(i,15)*b(15,j) &
654  + a(i,16)*b(16,j) &
655  + a(i,17)*b(17,j) &
656  + a(i,18)*b(18,j) &
657  + a(i,19)*b(19,j) &
658  + a(i,20)*b(20,j) &
659  + a(i,21)*b(21,j) &
660  + a(i,22)*b(22,j) &
661  + a(i,23)*b(23,j) &
662  + a(i,24)*b(24,j)
663  enddo
664  enddo
665  return
666  end subroutine mxf24
667 !-----------------------------------------------------------------------
668  subroutine mxm44_0(a, m, b, k, c, n)
669 
670 ! matrix multiply with a 4x4 pencil
671 
672  use kinds, only : dp
673  real(DP) :: a(m,k), b(k,n), c(m,n)
674  real(DP) :: s11, s12, s13, s14, s21, s22, s23, s24
675  real(DP) :: s31, s32, s33, s34, s41, s42, s43, s44
676 
677  mresid = iand(m,3)
678  nresid = iand(n,3)
679  m1 = m - mresid + 1
680  n1 = n - nresid + 1
681 
682  do i=1,m-mresid,4
683  do j=1,n-nresid,4
684  s11 = 0.0d0
685  s21 = 0.0d0
686  s31 = 0.0d0
687  s41 = 0.0d0
688  s12 = 0.0d0
689  s22 = 0.0d0
690  s32 = 0.0d0
691  s42 = 0.0d0
692  s13 = 0.0d0
693  s23 = 0.0d0
694  s33 = 0.0d0
695  s43 = 0.0d0
696  s14 = 0.0d0
697  s24 = 0.0d0
698  s34 = 0.0d0
699  s44 = 0.0d0
700  do l=1,k
701  s11 = s11 + a(i,l)*b(l,j)
702  s12 = s12 + a(i,l)*b(l,j+1)
703  s13 = s13 + a(i,l)*b(l,j+2)
704  s14 = s14 + a(i,l)*b(l,j+3)
705 
706  s21 = s21 + a(i+1,l)*b(l,j)
707  s22 = s22 + a(i+1,l)*b(l,j+1)
708  s23 = s23 + a(i+1,l)*b(l,j+2)
709  s24 = s24 + a(i+1,l)*b(l,j+3)
710 
711  s31 = s31 + a(i+2,l)*b(l,j)
712  s32 = s32 + a(i+2,l)*b(l,j+1)
713  s33 = s33 + a(i+2,l)*b(l,j+2)
714  s34 = s34 + a(i+2,l)*b(l,j+3)
715 
716  s41 = s41 + a(i+3,l)*b(l,j)
717  s42 = s42 + a(i+3,l)*b(l,j+1)
718  s43 = s43 + a(i+3,l)*b(l,j+2)
719  s44 = s44 + a(i+3,l)*b(l,j+3)
720  enddo
721  c(i,j) = s11
722  c(i,j+1) = s12
723  c(i,j+2) = s13
724  c(i,j+3) = s14
725 
726  c(i+1,j) = s21
727  c(i+2,j) = s31
728  c(i+3,j) = s41
729 
730  c(i+1,j+1) = s22
731  c(i+2,j+1) = s32
732  c(i+3,j+1) = s42
733 
734  c(i+1,j+2) = s23
735  c(i+2,j+2) = s33
736  c(i+3,j+2) = s43
737 
738  c(i+1,j+3) = s24
739  c(i+2,j+3) = s34
740  c(i+3,j+3) = s44
741  enddo
742  ! Residual when n is not multiple of 4
743  if (nresid /= 0) then
744  if (nresid == 1) then
745  s11 = 0.0d0
746  s21 = 0.0d0
747  s31 = 0.0d0
748  s41 = 0.0d0
749  do l=1,k
750  s11 = s11 + a(i,l)*b(l,n)
751  s21 = s21 + a(i+1,l)*b(l,n)
752  s31 = s31 + a(i+2,l)*b(l,n)
753  s41 = s41 + a(i+3,l)*b(l,n)
754  enddo
755  c(i,n) = s11
756  c(i+1,n) = s21
757  c(i+2,n) = s31
758  c(i+3,n) = s41
759  elseif (nresid == 2) then
760  s11 = 0.0d0
761  s21 = 0.0d0
762  s31 = 0.0d0
763  s41 = 0.0d0
764  s12 = 0.0d0
765  s22 = 0.0d0
766  s32 = 0.0d0
767  s42 = 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 
772  s21 = s21 + a(i+1,l)*b(l,j)
773  s22 = s22 + a(i+1,l)*b(l,j+1)
774 
775  s31 = s31 + a(i+2,l)*b(l,j)
776  s32 = s32 + a(i+2,l)*b(l,j+1)
777 
778  s41 = s41 + a(i+3,l)*b(l,j)
779  s42 = s42 + a(i+3,l)*b(l,j+1)
780  enddo
781  c(i,j) = s11
782  c(i,j+1) = s12
783 
784  c(i+1,j) = s21
785  c(i+2,j) = s31
786  c(i+3,j) = s41
787 
788  c(i+1,j+1) = s22
789  c(i+2,j+1) = s32
790  c(i+3,j+1) = s42
791  else
792  s11 = 0.0d0
793  s21 = 0.0d0
794  s31 = 0.0d0
795  s41 = 0.0d0
796  s12 = 0.0d0
797  s22 = 0.0d0
798  s32 = 0.0d0
799  s42 = 0.0d0
800  s13 = 0.0d0
801  s23 = 0.0d0
802  s33 = 0.0d0
803  s43 = 0.0d0
804  do l=1,k
805  s11 = s11 + a(i,l)*b(l,j)
806  s12 = s12 + a(i,l)*b(l,j+1)
807  s13 = s13 + a(i,l)*b(l,j+2)
808 
809  s21 = s21 + a(i+1,l)*b(l,j)
810  s22 = s22 + a(i+1,l)*b(l,j+1)
811  s23 = s23 + a(i+1,l)*b(l,j+2)
812 
813  s31 = s31 + a(i+2,l)*b(l,j)
814  s32 = s32 + a(i+2,l)*b(l,j+1)
815  s33 = s33 + a(i+2,l)*b(l,j+2)
816 
817  s41 = s41 + a(i+3,l)*b(l,j)
818  s42 = s42 + a(i+3,l)*b(l,j+1)
819  s43 = s43 + a(i+3,l)*b(l,j+2)
820  enddo
821  c(i,j) = s11
822  c(i+1,j) = s21
823  c(i+2,j) = s31
824  c(i+3,j) = s41
825  c(i,j+1) = s12
826  c(i+1,j+1) = s22
827  c(i+2,j+1) = s32
828  c(i+3,j+1) = s42
829  c(i,j+2) = s13
830  c(i+1,j+2) = s23
831  c(i+2,j+2) = s33
832  c(i+3,j+2) = s43
833  endif
834  endif
835  enddo
836 
837 ! Residual when m is not multiple of 4
838  if (mresid == 0) then
839  return
840  elseif (mresid == 1) then
841  do j=1,n-nresid,4
842  s11 = 0.0d0
843  s12 = 0.0d0
844  s13 = 0.0d0
845  s14 = 0.0d0
846  do l=1,k
847  s11 = s11 + a(m,l)*b(l,j)
848  s12 = s12 + a(m,l)*b(l,j+1)
849  s13 = s13 + a(m,l)*b(l,j+2)
850  s14 = s14 + a(m,l)*b(l,j+3)
851  enddo
852  c(m,j) = s11
853  c(m,j+1) = s12
854  c(m,j+2) = s13
855  c(m,j+3) = s14
856  enddo
857  ! mresid is 1, check nresid
858  if (nresid == 0) then
859  return
860  elseif (nresid == 1) then
861  s11 = 0.0d0
862  do l=1,k
863  s11 = s11 + a(m,l)*b(l,n)
864  enddo
865  c(m,n) = s11
866  return
867  elseif (nresid == 2) then
868  s11 = 0.0d0
869  s12 = 0.0d0
870  do l=1,k
871  s11 = s11 + a(m,l)*b(l,n-1)
872  s12 = s12 + a(m,l)*b(l,n)
873  enddo
874  c(m,n-1) = s11
875  c(m,n) = s12
876  return
877  else
878  s11 = 0.0d0
879  s12 = 0.0d0
880  s13 = 0.0d0
881  do l=1,k
882  s11 = s11 + a(m,l)*b(l,n-2)
883  s12 = s12 + a(m,l)*b(l,n-1)
884  s13 = s13 + a(m,l)*b(l,n)
885  enddo
886  c(m,n-2) = s11
887  c(m,n-1) = s12
888  c(m,n) = s13
889  return
890  endif
891  elseif (mresid == 2) then
892  do j=1,n-nresid,4
893  s11 = 0.0d0
894  s12 = 0.0d0
895  s13 = 0.0d0
896  s14 = 0.0d0
897  s21 = 0.0d0
898  s22 = 0.0d0
899  s23 = 0.0d0
900  s24 = 0.0d0
901  do l=1,k
902  s11 = s11 + a(m-1,l)*b(l,j)
903  s12 = s12 + a(m-1,l)*b(l,j+1)
904  s13 = s13 + a(m-1,l)*b(l,j+2)
905  s14 = s14 + a(m-1,l)*b(l,j+3)
906 
907  s21 = s21 + a(m,l)*b(l,j)
908  s22 = s22 + a(m,l)*b(l,j+1)
909  s23 = s23 + a(m,l)*b(l,j+2)
910  s24 = s24 + a(m,l)*b(l,j+3)
911  enddo
912  c(m-1,j) = s11
913  c(m-1,j+1) = s12
914  c(m-1,j+2) = s13
915  c(m-1,j+3) = s14
916  c(m,j) = s21
917  c(m,j+1) = s22
918  c(m,j+2) = s23
919  c(m,j+3) = s24
920  enddo
921  ! mresid is 2, check nresid
922  if (nresid == 0) then
923  return
924  elseif (nresid == 1) then
925  s11 = 0.0d0
926  s21 = 0.0d0
927  do l=1,k
928  s11 = s11 + a(m-1,l)*b(l,n)
929  s21 = s21 + a(m,l)*b(l,n)
930  enddo
931  c(m-1,n) = s11
932  c(m,n) = s21
933  return
934  elseif (nresid == 2) then
935  s11 = 0.0d0
936  s21 = 0.0d0
937  s12 = 0.0d0
938  s22 = 0.0d0
939  do l=1,k
940  s11 = s11 + a(m-1,l)*b(l,n-1)
941  s12 = s12 + a(m-1,l)*b(l,n)
942  s21 = s21 + a(m,l)*b(l,n-1)
943  s22 = s22 + a(m,l)*b(l,n)
944  enddo
945  c(m-1,n-1) = s11
946  c(m-1,n) = s12
947  c(m,n-1) = s21
948  c(m,n) = s22
949  return
950  else
951  s11 = 0.0d0
952  s21 = 0.0d0
953  s12 = 0.0d0
954  s22 = 0.0d0
955  s13 = 0.0d0
956  s23 = 0.0d0
957  do l=1,k
958  s11 = s11 + a(m-1,l)*b(l,n-2)
959  s12 = s12 + a(m-1,l)*b(l,n-1)
960  s13 = s13 + a(m-1,l)*b(l,n)
961  s21 = s21 + a(m,l)*b(l,n-2)
962  s22 = s22 + a(m,l)*b(l,n-1)
963  s23 = s23 + a(m,l)*b(l,n)
964  enddo
965  c(m-1,n-2) = s11
966  c(m-1,n-1) = s12
967  c(m-1,n) = s13
968  c(m,n-2) = s21
969  c(m,n-1) = s22
970  c(m,n) = s23
971  return
972  endif
973  else
974  ! mresid is 3
975  do j=1,n-nresid,4
976  s11 = 0.0d0
977  s21 = 0.0d0
978  s31 = 0.0d0
979 
980  s12 = 0.0d0
981  s22 = 0.0d0
982  s32 = 0.0d0
983 
984  s13 = 0.0d0
985  s23 = 0.0d0
986  s33 = 0.0d0
987 
988  s14 = 0.0d0
989  s24 = 0.0d0
990  s34 = 0.0d0
991 
992  do l=1,k
993  s11 = s11 + a(m-2,l)*b(l,j)
994  s12 = s12 + a(m-2,l)*b(l,j+1)
995  s13 = s13 + a(m-2,l)*b(l,j+2)
996  s14 = s14 + a(m-2,l)*b(l,j+3)
997 
998  s21 = s21 + a(m-1,l)*b(l,j)
999  s22 = s22 + a(m-1,l)*b(l,j+1)
1000  s23 = s23 + a(m-1,l)*b(l,j+2)
1001  s24 = s24 + a(m-1,l)*b(l,j+3)
1002 
1003  s31 = s31 + a(m,l)*b(l,j)
1004  s32 = s32 + a(m,l)*b(l,j+1)
1005  s33 = s33 + a(m,l)*b(l,j+2)
1006  s34 = s34 + a(m,l)*b(l,j+3)
1007  enddo
1008  c(m-2,j) = s11
1009  c(m-2,j+1) = s12
1010  c(m-2,j+2) = s13
1011  c(m-2,j+3) = s14
1012 
1013  c(m-1,j) = s21
1014  c(m-1,j+1) = s22
1015  c(m-1,j+2) = s23
1016  c(m-1,j+3) = s24
1017 
1018  c(m,j) = s31
1019  c(m,j+1) = s32
1020  c(m,j+2) = s33
1021  c(m,j+3) = s34
1022  enddo
1023  ! mresid is 3, check nresid
1024  if (nresid == 0) then
1025  return
1026  elseif (nresid == 1) then
1027  s11 = 0.0d0
1028  s21 = 0.0d0
1029  s31 = 0.0d0
1030  do l=1,k
1031  s11 = s11 + a(m-2,l)*b(l,n)
1032  s21 = s21 + a(m-1,l)*b(l,n)
1033  s31 = s31 + a(m,l)*b(l,n)
1034  enddo
1035  c(m-2,n) = s11
1036  c(m-1,n) = s21
1037  c(m,n) = s31
1038  return
1039  elseif (nresid == 2) then
1040  s11 = 0.0d0
1041  s21 = 0.0d0
1042  s31 = 0.0d0
1043  s12 = 0.0d0
1044  s22 = 0.0d0
1045  s32 = 0.0d0
1046  do l=1,k
1047  s11 = s11 + a(m-2,l)*b(l,n-1)
1048  s12 = s12 + a(m-2,l)*b(l,n)
1049  s21 = s21 + a(m-1,l)*b(l,n-1)
1050  s22 = s22 + a(m-1,l)*b(l,n)
1051  s31 = s31 + a(m,l)*b(l,n-1)
1052  s32 = s32 + a(m,l)*b(l,n)
1053  enddo
1054  c(m-2,n-1) = s11
1055  c(m-2,n) = s12
1056  c(m-1,n-1) = s21
1057  c(m-1,n) = s22
1058  c(m,n-1) = s31
1059  c(m,n) = s32
1060  return
1061  else
1062  s11 = 0.0d0
1063  s21 = 0.0d0
1064  s31 = 0.0d0
1065  s12 = 0.0d0
1066  s22 = 0.0d0
1067  s32 = 0.0d0
1068  s13 = 0.0d0
1069  s23 = 0.0d0
1070  s33 = 0.0d0
1071  do l=1,k
1072  s11 = s11 + a(m-2,l)*b(l,n-2)
1073  s12 = s12 + a(m-2,l)*b(l,n-1)
1074  s13 = s13 + a(m-2,l)*b(l,n)
1075  s21 = s21 + a(m-1,l)*b(l,n-2)
1076  s22 = s22 + a(m-1,l)*b(l,n-1)
1077  s23 = s23 + a(m-1,l)*b(l,n)
1078  s31 = s31 + a(m,l)*b(l,n-2)
1079  s32 = s32 + a(m,l)*b(l,n-1)
1080  s33 = s33 + a(m,l)*b(l,n)
1081  enddo
1082  c(m-2,n-2) = s11
1083  c(m-2,n-1) = s12
1084  c(m-2,n) = s13
1085  c(m-1,n-2) = s21
1086  c(m-1,n-1) = s22
1087  c(m-1,n) = s23
1088  c(m,n-2) = s31
1089  c(m,n-1) = s32
1090  c(m,n) = s33
1091  return
1092  endif
1093  endif
1094 
1095  return
1096  end subroutine mxm44_0
1097 !-----------------------------------------------------------------------
1098  subroutine mxm44_2(a, m, b, k, c, n)
1099  use kinds, only : dp
1100  real(DP) :: a(m,2), b(2,n), c(m,n)
1101 
1102  nresid = iand(n,3)
1103  n1 = n - nresid + 1
1104 
1105  do j=1,n-nresid,4
1106  do i=1,m
1107  c(i,j) = a(i,1)*b(1,j) &
1108  + a(i,2)*b(2,j)
1109  c(i,j+1) = a(i,1)*b(1,j+1) &
1110  + a(i,2)*b(2,j+1)
1111  c(i,j+2) = a(i,1)*b(1,j+2) &
1112  + a(i,2)*b(2,j+2)
1113  c(i,j+3) = a(i,1)*b(1,j+3) &
1114  + a(i,2)*b(2,j+3)
1115  enddo
1116  enddo
1117  if (nresid == 0) then
1118  return
1119  elseif (nresid == 1) then
1120  do i=1,m
1121  c(i,n) = a(i,1)*b(1,n) &
1122  + a(i,2)*b(2,n)
1123  enddo
1124  elseif (nresid == 2) then
1125  do i=1,m
1126  c(i,n-1) = a(i,1)*b(1,n-1) &
1127  + a(i,2)*b(2,n-1)
1128  c(i,n) = a(i,1)*b(1,n) &
1129  + a(i,2)*b(2,n)
1130  enddo
1131  else
1132  do i=1,m
1133  c(i,n-2) = a(i,1)*b(1,n-2) &
1134  + a(i,2)*b(2,n-2)
1135  c(i,n-1) = a(i,1)*b(1,n-1) &
1136  + a(i,2)*b(2,n-1)
1137  c(i,n) = a(i,1)*b(1,n) &
1138  + a(i,2)*b(2,n)
1139  enddo
1140  endif
1141 
1142  return
1143  end subroutine mxm44_2
1144 !-----------------------------------------------------------------------
1145  subroutine initab(a,b,n)
1146  use kinds, only : dp
1147  real(DP) :: a(1),b(1)
1148  do i=1,n-1
1149  x = i
1150  k = mod(i,19) + 2
1151  l = mod(i,17) + 5
1152  m = mod(i,31) + 3
1153  a(i) = -.25*(a(i)+a(i+1)) + (x*x + k + l)/(x*x+m)
1154  b(i) = -.25*(b(i)+b(i+1)) + (x*x + k + m)/(x*x+l)
1155  enddo
1156  a(n) = -.25*(a(n)+a(n)) + (x*x + k + l)/(x*x+m)
1157  b(n) = -.25*(b(n)+b(n)) + (x*x + k + m)/(x*x+l)
1158  return
1159  end subroutine initab
1160 !-----------------------------------------------------------------------
1161  subroutine mxms(a,n1,b,n2,c,n3)
1162 !----------------------------------------------------------------------
1163 
1164 ! Matrix-vector product routine.
1165 ! NOTE: Use assembly coded routine if available.
1166 
1167 !---------------------------------------------------------------------
1168  use kinds, only : dp
1169  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
1170 
1171  n0=n1*n3
1172  DO 10 i=1,n0
1173  c(i,1)=0.
1174  10 END DO
1175  DO 100 j=1,n3
1176  DO 100 k=1,n2
1177  bb=b(k,j)
1178  DO 100 i=1,n1
1179  c(i,j)=c(i,j)+a(i,k)*bb
1180  100 END DO
1181  return
1182  end subroutine mxms
1183 !-----------------------------------------------------------------------
1184  subroutine mxmu4(a,n1,b,n2,c,n3)
1185 !----------------------------------------------------------------------
1186 
1187 ! Matrix-vector product routine.
1188 ! NOTE: Use assembly coded routine if available.
1189 
1190 !---------------------------------------------------------------------
1191  use kinds, only : dp
1192  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
1193 
1194  n0=n1*n3
1195  DO 10 i=1,n0
1196  c(i,1)=0.
1197  10 END DO
1198  i1 = n1 - mod(n1,4) + 1
1199  DO 100 j=1,n3
1200  DO 100 k=1,n2
1201  bb=b(k,j)
1202  DO i=1,n1-3,4
1203  c(i ,j)=c(i ,j)+a(i ,k)*bb
1204  c(i+1,j)=c(i+1,j)+a(i+1,k)*bb
1205  c(i+2,j)=c(i+2,j)+a(i+2,k)*bb
1206  c(i+3,j)=c(i+3,j)+a(i+3,k)*bb
1207  enddo
1208  DO i=i1,n1
1209  c(i ,j)=c(i ,j)+a(i ,k)*bb
1210  enddo
1211  100 END DO
1212  return
1213  end subroutine mxmu4
1214 !-----------------------------------------------------------------------
1215  subroutine madd (a,n1,b,n2,c,n3)
1216 
1217  use kinds, only : dp
1218  real(DP) :: a(n1,n2),b(n2,n3),c(n1,n3)
1219 
1220  do j=1,n3
1221  do i=1,n1
1222  c(i,j) = a(i,j)+b(i,j)
1223  enddo
1224  enddo
1225 
1226  return
1227  end subroutine madd
1228 !-----------------------------------------------------------------------
1229  subroutine mxmur2(a,n1,b,n2,c,n3)
1230 !----------------------------------------------------------------------
1231 
1232 ! Matrix-vector product routine.
1233 ! NOTE: Use assembly coded routine if available.
1234 
1235 !---------------------------------------------------------------------
1236  use kinds, only : dp
1237  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
1238 
1239  if (n2 <= 8) then
1240  if (n2 == 1) then
1241  call mxmur2_1(a,n1,b,n2,c,n3)
1242  elseif (n2 == 2) then
1243  call mxmur2_2(a,n1,b,n2,c,n3)
1244  elseif (n2 == 3) then
1245  call mxmur2_3(a,n1,b,n2,c,n3)
1246  elseif (n2 == 4) then
1247  call mxmur2_4(a,n1,b,n2,c,n3)
1248  elseif (n2 == 5) then
1249  call mxmur2_5(a,n1,b,n2,c,n3)
1250  elseif (n2 == 6) then
1251  call mxmur2_6(a,n1,b,n2,c,n3)
1252  elseif (n2 == 7) then
1253  call mxmur2_7(a,n1,b,n2,c,n3)
1254  else
1255  call mxmur2_8(a,n1,b,n2,c,n3)
1256  endif
1257  elseif (n2 <= 16) then
1258  if (n2 == 9) then
1259  call mxmur2_9(a,n1,b,n2,c,n3)
1260  elseif (n2 == 10) then
1261  call mxmur2_10(a,n1,b,n2,c,n3)
1262  elseif (n2 == 11) then
1263  call mxmur2_11(a,n1,b,n2,c,n3)
1264  elseif (n2 == 12) then
1265  call mxmur2_12(a,n1,b,n2,c,n3)
1266  elseif (n2 == 13) then
1267  call mxmur2_13(a,n1,b,n2,c,n3)
1268  elseif (n2 == 14) then
1269  call mxmur2_14(a,n1,b,n2,c,n3)
1270  elseif (n2 == 15) then
1271  call mxmur2_15(a,n1,b,n2,c,n3)
1272  else
1273  call mxmur2_16(a,n1,b,n2,c,n3)
1274  endif
1275  else
1276  n0=n1*n3
1277  DO 10 i=1,n0
1278  c(i,1)=0.
1279  10 END DO
1280  DO 100 j=1,n3
1281  DO 100 k=1,n2
1282  bb=b(k,j)
1283  DO 100 i=1,n1
1284  c(i,j)=c(i,j)+a(i,k)*bb
1285  100 END DO
1286  endif
1287  return
1288  end subroutine mxmur2
1289 
1290  subroutine mxmur2_1(a,n1,b,n2,c,n3)
1291 
1292  use kinds, only : dp
1293  real(DP) :: a(n1,1),b(1,n3),c(n1,n3)
1294 
1295  do j=1,n3
1296  do i=1,n1
1297  c(i,j) = a(i,1)*b(1,j)
1298  enddo
1299  enddo
1300  return
1301  end subroutine mxmur2_1
1302  subroutine mxmur2_2(a,n1,b,n2,c,n3)
1303 
1304  use kinds, only : dp
1305  real(DP) :: a(n1,2),b(2,n3),c(n1,n3)
1306 
1307  do j=1,n3
1308  do i=1,n1
1309  c(i,j) = a(i,1)*b(1,j) &
1310  + a(i,2)*b(2,j)
1311  enddo
1312  enddo
1313  return
1314  end subroutine mxmur2_2
1315  subroutine mxmur2_3(a,n1,b,n2,c,n3)
1316 
1317  use kinds, only : dp
1318  real(DP) :: a(n1,3),b(3,n3),c(n1,n3)
1319 
1320  do j=1,n3
1321  do i=1,n1
1322  c(i,j) = a(i,1)*b(1,j) &
1323  + a(i,2)*b(2,j) &
1324  + a(i,3)*b(3,j)
1325  enddo
1326  enddo
1327  return
1328  end subroutine mxmur2_3
1329  subroutine mxmur2_4(a,n1,b,n2,c,n3)
1330 
1331  use kinds, only : dp
1332  real(DP) :: a(n1,4),b(4,n3),c(n1,n3)
1333 
1334  do j=1,n3
1335  do i=1,n1
1336  c(i,j) = a(i,1)*b(1,j) &
1337  + a(i,2)*b(2,j) &
1338  + a(i,3)*b(3,j) &
1339  + a(i,4)*b(4,j)
1340  enddo
1341  enddo
1342  return
1343  end subroutine mxmur2_4
1344  subroutine mxmur2_5(a,n1,b,n2,c,n3)
1345 
1346  use kinds, only : dp
1347  real(DP) :: a(n1,5),b(5,n3),c(n1,n3)
1348 
1349  do j=1,n3
1350  do i=1,n1
1351  c(i,j) = a(i,1)*b(1,j) &
1352  + a(i,2)*b(2,j) &
1353  + a(i,3)*b(3,j) &
1354  + a(i,4)*b(4,j) &
1355  + a(i,5)*b(5,j)
1356  enddo
1357  enddo
1358  return
1359  end subroutine mxmur2_5
1360  subroutine mxmur2_6(a,n1,b,n2,c,n3)
1361 
1362  use kinds, only : dp
1363  real(DP) :: a(n1,6),b(6,n3),c(n1,n3)
1364 
1365  do j=1,n3
1366  do i=1,n1
1367  c(i,j) = a(i,1)*b(1,j) &
1368  + a(i,2)*b(2,j) &
1369  + a(i,3)*b(3,j) &
1370  + a(i,4)*b(4,j) &
1371  + a(i,5)*b(5,j) &
1372  + a(i,6)*b(6,j)
1373  enddo
1374  enddo
1375  return
1376  end subroutine mxmur2_6
1377  subroutine mxmur2_7(a,n1,b,n2,c,n3)
1378 
1379  use kinds, only : dp
1380  real(DP) :: a(n1,7),b(7,n3),c(n1,n3)
1381 
1382  do j=1,n3
1383  do i=1,n1
1384  c(i,j) = a(i,1)*b(1,j) &
1385  + a(i,2)*b(2,j) &
1386  + a(i,3)*b(3,j) &
1387  + a(i,4)*b(4,j) &
1388  + a(i,5)*b(5,j) &
1389  + a(i,6)*b(6,j) &
1390  + a(i,7)*b(7,j)
1391  enddo
1392  enddo
1393  return
1394  end subroutine mxmur2_7
1395  subroutine mxmur2_8(a,n1,b,n2,c,n3)
1396 
1397  use kinds, only : dp
1398  real(DP) :: a(n1,8),b(8,n3),c(n1,n3)
1399 
1400  do j=1,n3
1401  do i=1,n1
1402  c(i,j) = a(i,1)*b(1,j) &
1403  + a(i,2)*b(2,j) &
1404  + a(i,3)*b(3,j) &
1405  + a(i,4)*b(4,j) &
1406  + a(i,5)*b(5,j) &
1407  + a(i,6)*b(6,j) &
1408  + a(i,7)*b(7,j) &
1409  + a(i,8)*b(8,j)
1410  enddo
1411  enddo
1412  return
1413  end subroutine mxmur2_8
1414  subroutine mxmur2_9(a,n1,b,n2,c,n3)
1415 
1416  use kinds, only : dp
1417  real(DP) :: a(n1,9),b(9,n3),c(n1,n3)
1418 
1419  do j=1,n3
1420  do i=1,n1
1421  c(i,j) = a(i,1)*b(1,j) &
1422  + a(i,2)*b(2,j) &
1423  + a(i,3)*b(3,j) &
1424  + a(i,4)*b(4,j) &
1425  + a(i,5)*b(5,j) &
1426  + a(i,6)*b(6,j) &
1427  + a(i,7)*b(7,j) &
1428  + a(i,8)*b(8,j) &
1429  + a(i,9)*b(9,j)
1430  enddo
1431  enddo
1432  return
1433  end subroutine mxmur2_9
1434  subroutine mxmur2_10(a,n1,b,n2,c,n3)
1435 
1436  use kinds, only : dp
1437  real(DP) :: a(n1,10),b(10,n3),c(n1,n3)
1438 
1439  do j=1,n3
1440  do i=1,n1
1441  c(i,j) = a(i,1)*b(1,j) &
1442  + a(i,2)*b(2,j) &
1443  + a(i,3)*b(3,j) &
1444  + a(i,4)*b(4,j) &
1445  + a(i,5)*b(5,j) &
1446  + a(i,6)*b(6,j) &
1447  + a(i,7)*b(7,j) &
1448  + a(i,8)*b(8,j) &
1449  + a(i,9)*b(9,j) &
1450  + a(i,10)*b(10,j)
1451  enddo
1452  enddo
1453  return
1454  end subroutine mxmur2_10
1455  subroutine mxmur2_11(a,n1,b,n2,c,n3)
1456 
1457  use kinds, only : dp
1458  real(DP) :: a(n1,11),b(11,n3),c(n1,n3)
1459 
1460  do j=1,n3
1461  do i=1,n1
1462  c(i,j) = a(i,1)*b(1,j) &
1463  + a(i,2)*b(2,j) &
1464  + a(i,3)*b(3,j) &
1465  + a(i,4)*b(4,j) &
1466  + a(i,5)*b(5,j) &
1467  + a(i,6)*b(6,j) &
1468  + a(i,7)*b(7,j) &
1469  + a(i,8)*b(8,j) &
1470  + a(i,9)*b(9,j) &
1471  + a(i,10)*b(10,j) &
1472  + a(i,11)*b(11,j)
1473  enddo
1474  enddo
1475  return
1476  end subroutine mxmur2_11
1477  subroutine mxmur2_12(a,n1,b,n2,c,n3)
1478 
1479  use kinds, only : dp
1480  real(DP) :: a(n1,12),b(12,n3),c(n1,n3)
1481 
1482  do j=1,n3
1483  do i=1,n1
1484  c(i,j) = a(i,1)*b(1,j) &
1485  + a(i,2)*b(2,j) &
1486  + a(i,3)*b(3,j) &
1487  + a(i,4)*b(4,j) &
1488  + a(i,5)*b(5,j) &
1489  + a(i,6)*b(6,j) &
1490  + a(i,7)*b(7,j) &
1491  + a(i,8)*b(8,j) &
1492  + a(i,9)*b(9,j) &
1493  + a(i,10)*b(10,j) &
1494  + a(i,11)*b(11,j) &
1495  + a(i,12)*b(12,j)
1496  enddo
1497  enddo
1498  return
1499  end subroutine mxmur2_12
1500  subroutine mxmur2_13(a,n1,b,n2,c,n3)
1501 
1502  use kinds, only : dp
1503  real(DP) :: a(n1,13),b(13,n3),c(n1,n3)
1504 
1505  do j=1,n3
1506  do i=1,n1
1507  c(i,j) = a(i,1)*b(1,j) &
1508  + a(i,2)*b(2,j) &
1509  + a(i,3)*b(3,j) &
1510  + a(i,4)*b(4,j) &
1511  + a(i,5)*b(5,j) &
1512  + a(i,6)*b(6,j) &
1513  + a(i,7)*b(7,j) &
1514  + a(i,8)*b(8,j) &
1515  + a(i,9)*b(9,j) &
1516  + a(i,10)*b(10,j) &
1517  + a(i,11)*b(11,j) &
1518  + a(i,12)*b(12,j) &
1519  + a(i,13)*b(13,j)
1520  enddo
1521  enddo
1522  return
1523  end subroutine mxmur2_13
1524  subroutine mxmur2_14(a,n1,b,n2,c,n3)
1525 
1526  use kinds, only : dp
1527  real(DP) :: a(n1,14),b(14,n3),c(n1,n3)
1528 
1529  do j=1,n3
1530  do i=1,n1
1531  c(i,j) = a(i,1)*b(1,j) &
1532  + a(i,2)*b(2,j) &
1533  + a(i,3)*b(3,j) &
1534  + a(i,4)*b(4,j) &
1535  + a(i,5)*b(5,j) &
1536  + a(i,6)*b(6,j) &
1537  + a(i,7)*b(7,j) &
1538  + a(i,8)*b(8,j) &
1539  + a(i,9)*b(9,j) &
1540  + a(i,10)*b(10,j) &
1541  + a(i,11)*b(11,j) &
1542  + a(i,12)*b(12,j) &
1543  + a(i,13)*b(13,j) &
1544  + a(i,14)*b(14,j)
1545  enddo
1546  enddo
1547  return
1548  end subroutine mxmur2_14
1549  subroutine mxmur2_15(a,n1,b,n2,c,n3)
1550 
1551  use kinds, only : dp
1552  real(DP) :: a(n1,15),b(15,n3),c(n1,n3)
1553 
1554  do j=1,n3
1555  do i=1,n1
1556  c(i,j) = a(i,1)*b(1,j) &
1557  + a(i,2)*b(2,j) &
1558  + a(i,3)*b(3,j) &
1559  + a(i,4)*b(4,j) &
1560  + a(i,5)*b(5,j) &
1561  + a(i,6)*b(6,j) &
1562  + a(i,7)*b(7,j) &
1563  + a(i,8)*b(8,j) &
1564  + a(i,9)*b(9,j) &
1565  + a(i,10)*b(10,j) &
1566  + a(i,11)*b(11,j) &
1567  + a(i,12)*b(12,j) &
1568  + a(i,13)*b(13,j) &
1569  + a(i,14)*b(14,j) &
1570  + a(i,15)*b(15,j)
1571  enddo
1572  enddo
1573  return
1574  end subroutine mxmur2_15
1575  subroutine mxmur2_16(a,n1,b,n2,c,n3)
1576 
1577  use kinds, only : dp
1578  real(DP) :: a(n1,16),b(16,n3),c(n1,n3)
1579 
1580  do j=1,n3
1581  do i=1,n1
1582  c(i,j) = a(i,1)*b(1,j) &
1583  + a(i,2)*b(2,j) &
1584  + a(i,3)*b(3,j) &
1585  + a(i,4)*b(4,j) &
1586  + a(i,5)*b(5,j) &
1587  + a(i,6)*b(6,j) &
1588  + a(i,7)*b(7,j) &
1589  + a(i,8)*b(8,j) &
1590  + a(i,9)*b(9,j) &
1591  + a(i,10)*b(10,j) &
1592  + a(i,11)*b(11,j) &
1593  + a(i,12)*b(12,j) &
1594  + a(i,13)*b(13,j) &
1595  + a(i,14)*b(14,j) &
1596  + a(i,15)*b(15,j) &
1597  + a(i,16)*b(16,j)
1598  enddo
1599  enddo
1600  return
1601  end subroutine mxmur2_16
1602 !-----------------------------------------------------------------------
1603  subroutine mxmur3(a,n1,b,n2,c,n3)
1604 !----------------------------------------------------------------------
1605 
1606 ! Matrix-vector product routine.
1607 ! NOTE: Use assembly coded routine if available.
1608 
1609 !---------------------------------------------------------------------
1610  use kinds, only : dp
1611  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
1612 
1613  n0=n1*n3
1614  DO 10 i=1,n0
1615  c(i,1)=0.
1616  10 END DO
1617  if (n3 <= 8) then
1618  if (n3 == 1) then
1619  call mxmur3_1(a,n1,b,n2,c,n3)
1620  elseif (n3 == 2) then
1621  call mxmur3_2(a,n1,b,n2,c,n3)
1622  elseif (n3 == 3) then
1623  call mxmur3_3(a,n1,b,n2,c,n3)
1624  elseif (n3 == 4) then
1625  call mxmur3_4(a,n1,b,n2,c,n3)
1626  elseif (n3 == 5) then
1627  call mxmur3_5(a,n1,b,n2,c,n3)
1628  elseif (n3 == 6) then
1629  call mxmur3_6(a,n1,b,n2,c,n3)
1630  elseif (n3 == 7) then
1631  call mxmur3_7(a,n1,b,n2,c,n3)
1632  else
1633  call mxmur3_8(a,n1,b,n2,c,n3)
1634  endif
1635  elseif (n3 <= 16) then
1636  if (n3 == 9) then
1637  call mxmur3_9(a,n1,b,n2,c,n3)
1638  elseif (n3 == 10) then
1639  call mxmur3_10(a,n1,b,n2,c,n3)
1640  elseif (n3 == 11) then
1641  call mxmur3_11(a,n1,b,n2,c,n3)
1642  elseif (n3 == 12) then
1643  call mxmur3_12(a,n1,b,n2,c,n3)
1644  elseif (n3 == 13) then
1645  call mxmur3_13(a,n1,b,n2,c,n3)
1646  elseif (n3 == 14) then
1647  call mxmur3_14(a,n1,b,n2,c,n3)
1648  elseif (n3 == 15) then
1649  call mxmur3_15(a,n1,b,n2,c,n3)
1650  else
1651  call mxmur3_16(a,n1,b,n2,c,n3)
1652  endif
1653  else
1654  DO 100 j=1,n3
1655  DO 100 k=1,n2
1656  bb=b(k,j)
1657  DO 100 i=1,n1
1658  c(i,j)=c(i,j)+a(i,k)*bb
1659  100 END DO
1660  endif
1661  return
1662  end subroutine mxmur3
1663 
1664  subroutine mxmur3_16(a,n1,b,n2,c,n3)
1665  use kinds, only : dp
1666  real(DP) :: a(n1,n2),b(n2,16),c(n1,16)
1667 
1668  do k=1,n2
1669  tmp1 = b(k, 1)
1670  tmp2 = b(k, 2)
1671  tmp3 = b(k, 3)
1672  tmp4 = b(k, 4)
1673  tmp5 = b(k, 5)
1674  tmp6 = b(k, 6)
1675  tmp7 = b(k, 7)
1676  tmp8 = b(k, 8)
1677  tmp9 = b(k, 9)
1678  tmp10 = b(k,10)
1679  tmp11 = b(k,11)
1680  tmp12 = b(k,12)
1681  tmp13 = b(k,13)
1682  tmp14 = b(k,14)
1683  tmp15 = b(k,15)
1684  tmp16 = b(k,16)
1685  do i=1,n1
1686  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1687  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1688  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1689  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1690  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1691  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1692  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1693  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1694  c(i, 9) = c(i, 9) + a(i,k) * tmp9
1695  c(i,10) = c(i,10) + a(i,k) * tmp10
1696  c(i,11) = c(i,11) + a(i,k) * tmp11
1697  c(i,12) = c(i,12) + a(i,k) * tmp12
1698  c(i,13) = c(i,13) + a(i,k) * tmp13
1699  c(i,14) = c(i,14) + a(i,k) * tmp14
1700  c(i,15) = c(i,15) + a(i,k) * tmp15
1701  c(i,16) = c(i,16) + a(i,k) * tmp16
1702  enddo
1703 
1704  enddo
1705 
1706  return
1707  end subroutine mxmur3_16
1708  subroutine mxmur3_15(a,n1,b,n2,c,n3)
1709  use kinds, only : dp
1710  real(DP) :: a(n1,n2),b(n2,15),c(n1,15)
1711 
1712  do k=1,n2
1713  tmp1 = b(k, 1)
1714  tmp2 = b(k, 2)
1715  tmp3 = b(k, 3)
1716  tmp4 = b(k, 4)
1717  tmp5 = b(k, 5)
1718  tmp6 = b(k, 6)
1719  tmp7 = b(k, 7)
1720  tmp8 = b(k, 8)
1721  tmp9 = b(k, 9)
1722  tmp10 = b(k,10)
1723  tmp11 = b(k,11)
1724  tmp12 = b(k,12)
1725  tmp13 = b(k,13)
1726  tmp14 = b(k,14)
1727  tmp15 = b(k,15)
1728  do i=1,n1
1729  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1730  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1731  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1732  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1733  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1734  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1735  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1736  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1737  c(i, 9) = c(i, 9) + a(i,k) * tmp9
1738  c(i,10) = c(i,10) + a(i,k) * tmp10
1739  c(i,11) = c(i,11) + a(i,k) * tmp11
1740  c(i,12) = c(i,12) + a(i,k) * tmp12
1741  c(i,13) = c(i,13) + a(i,k) * tmp13
1742  c(i,14) = c(i,14) + a(i,k) * tmp14
1743  c(i,15) = c(i,15) + a(i,k) * tmp15
1744  enddo
1745 
1746  enddo
1747 
1748  return
1749  end subroutine mxmur3_15
1750  subroutine mxmur3_14(a,n1,b,n2,c,n3)
1751  use kinds, only : dp
1752  real(DP) :: a(n1,n2),b(n2,14),c(n1,14)
1753 
1754  do k=1,n2
1755  tmp1 = b(k, 1)
1756  tmp2 = b(k, 2)
1757  tmp3 = b(k, 3)
1758  tmp4 = b(k, 4)
1759  tmp5 = b(k, 5)
1760  tmp6 = b(k, 6)
1761  tmp7 = b(k, 7)
1762  tmp8 = b(k, 8)
1763  tmp9 = b(k, 9)
1764  tmp10 = b(k,10)
1765  tmp11 = b(k,11)
1766  tmp12 = b(k,12)
1767  tmp13 = b(k,13)
1768  tmp14 = b(k,14)
1769  do i=1,n1
1770  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1771  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1772  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1773  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1774  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1775  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1776  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1777  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1778  c(i, 9) = c(i, 9) + a(i,k) * tmp9
1779  c(i,10) = c(i,10) + a(i,k) * tmp10
1780  c(i,11) = c(i,11) + a(i,k) * tmp11
1781  c(i,12) = c(i,12) + a(i,k) * tmp12
1782  c(i,13) = c(i,13) + a(i,k) * tmp13
1783  c(i,14) = c(i,14) + a(i,k) * tmp14
1784  enddo
1785 
1786  enddo
1787 
1788  return
1789  end subroutine mxmur3_14
1790  subroutine mxmur3_13(a,n1,b,n2,c,n3)
1791  use kinds, only : dp
1792  real(DP) :: a(n1,n2),b(n2,13),c(n1,13)
1793 
1794  do k=1,n2
1795  tmp1 = b(k, 1)
1796  tmp2 = b(k, 2)
1797  tmp3 = b(k, 3)
1798  tmp4 = b(k, 4)
1799  tmp5 = b(k, 5)
1800  tmp6 = b(k, 6)
1801  tmp7 = b(k, 7)
1802  tmp8 = b(k, 8)
1803  tmp9 = b(k, 9)
1804  tmp10 = b(k,10)
1805  tmp11 = b(k,11)
1806  tmp12 = b(k,12)
1807  tmp13 = b(k,13)
1808  do i=1,n1
1809  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1810  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1811  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1812  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1813  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1814  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1815  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1816  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1817  c(i, 9) = c(i, 9) + a(i,k) * tmp9
1818  c(i,10) = c(i,10) + a(i,k) * tmp10
1819  c(i,11) = c(i,11) + a(i,k) * tmp11
1820  c(i,12) = c(i,12) + a(i,k) * tmp12
1821  c(i,13) = c(i,13) + a(i,k) * tmp13
1822  enddo
1823 
1824  enddo
1825 
1826  return
1827  end subroutine mxmur3_13
1828  subroutine mxmur3_12(a,n1,b,n2,c,n3)
1829  use kinds, only : dp
1830  real(DP) :: a(n1,n2),b(n2,12),c(n1,12)
1831 
1832  do k=1,n2
1833  tmp1 = b(k, 1)
1834  tmp2 = b(k, 2)
1835  tmp3 = b(k, 3)
1836  tmp4 = b(k, 4)
1837  tmp5 = b(k, 5)
1838  tmp6 = b(k, 6)
1839  tmp7 = b(k, 7)
1840  tmp8 = b(k, 8)
1841  tmp9 = b(k, 9)
1842  tmp10 = b(k,10)
1843  tmp11 = b(k,11)
1844  tmp12 = b(k,12)
1845  do i=1,n1
1846  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1847  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1848  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1849  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1850  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1851  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1852  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1853  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1854  c(i, 9) = c(i, 9) + a(i,k) * tmp9
1855  c(i,10) = c(i,10) + a(i,k) * tmp10
1856  c(i,11) = c(i,11) + a(i,k) * tmp11
1857  c(i,12) = c(i,12) + a(i,k) * tmp12
1858  enddo
1859 
1860  enddo
1861 
1862  return
1863  end subroutine mxmur3_12
1864  subroutine mxmur3_11(a,n1,b,n2,c,n3)
1865  use kinds, only : dp
1866  real(DP) :: a(n1,n2),b(n2,11),c(n1,11)
1867 
1868  do k=1,n2
1869  tmp1 = b(k, 1)
1870  tmp2 = b(k, 2)
1871  tmp3 = b(k, 3)
1872  tmp4 = b(k, 4)
1873  tmp5 = b(k, 5)
1874  tmp6 = b(k, 6)
1875  tmp7 = b(k, 7)
1876  tmp8 = b(k, 8)
1877  tmp9 = b(k, 9)
1878  tmp10 = b(k,10)
1879  tmp11 = b(k,11)
1880  do i=1,n1
1881  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1882  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1883  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1884  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1885  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1886  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1887  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1888  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1889  c(i, 9) = c(i, 9) + a(i,k) * tmp9
1890  c(i,10) = c(i,10) + a(i,k) * tmp10
1891  c(i,11) = c(i,11) + a(i,k) * tmp11
1892  enddo
1893 
1894  enddo
1895 
1896  return
1897  end subroutine mxmur3_11
1898  subroutine mxmur3_10(a,n1,b,n2,c,n3)
1899  use kinds, only : dp
1900  real(DP) :: a(n1,n2),b(n2,10),c(n1,10)
1901 
1902  do k=1,n2
1903  tmp1 = b(k, 1)
1904  tmp2 = b(k, 2)
1905  tmp3 = b(k, 3)
1906  tmp4 = b(k, 4)
1907  tmp5 = b(k, 5)
1908  tmp6 = b(k, 6)
1909  tmp7 = b(k, 7)
1910  tmp8 = b(k, 8)
1911  tmp9 = b(k, 9)
1912  tmp10 = b(k,10)
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  c(i,10) = c(i,10) + a(i,k) * tmp10
1924  enddo
1925 
1926  enddo
1927 
1928  return
1929  end subroutine mxmur3_10
1930  subroutine mxmur3_9(a,n1,b,n2,c,n3)
1931  use kinds, only : dp
1932  real(DP) :: a(n1,n2),b(n2,9),c(n1,9)
1933 
1934  do k=1,n2
1935  tmp1 = b(k, 1)
1936  tmp2 = b(k, 2)
1937  tmp3 = b(k, 3)
1938  tmp4 = b(k, 4)
1939  tmp5 = b(k, 5)
1940  tmp6 = b(k, 6)
1941  tmp7 = b(k, 7)
1942  tmp8 = b(k, 8)
1943  tmp9 = b(k, 9)
1944  do i=1,n1
1945  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1946  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1947  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1948  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1949  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1950  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1951  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1952  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1953  c(i, 9) = c(i, 9) + a(i,k) * tmp9
1954  enddo
1955 
1956  enddo
1957 
1958  return
1959  end subroutine mxmur3_9
1960  subroutine mxmur3_8(a,n1,b,n2,c,n3)
1961  use kinds, only : dp
1962  real(DP) :: a(n1,n2),b(n2,8),c(n1,8)
1963 
1964  do k=1,n2
1965  tmp1 = b(k, 1)
1966  tmp2 = b(k, 2)
1967  tmp3 = b(k, 3)
1968  tmp4 = b(k, 4)
1969  tmp5 = b(k, 5)
1970  tmp6 = b(k, 6)
1971  tmp7 = b(k, 7)
1972  tmp8 = b(k, 8)
1973  do i=1,n1
1974  c(i, 1) = c(i, 1) + a(i,k) * tmp1
1975  c(i, 2) = c(i, 2) + a(i,k) * tmp2
1976  c(i, 3) = c(i, 3) + a(i,k) * tmp3
1977  c(i, 4) = c(i, 4) + a(i,k) * tmp4
1978  c(i, 5) = c(i, 5) + a(i,k) * tmp5
1979  c(i, 6) = c(i, 6) + a(i,k) * tmp6
1980  c(i, 7) = c(i, 7) + a(i,k) * tmp7
1981  c(i, 8) = c(i, 8) + a(i,k) * tmp8
1982  enddo
1983 
1984  enddo
1985 
1986  return
1987  end subroutine mxmur3_8
1988  subroutine mxmur3_7(a,n1,b,n2,c,n3)
1989  use kinds, only : dp
1990  real(DP) :: a(n1,n2),b(n2,7),c(n1,7)
1991 
1992  do k=1,n2
1993  tmp1 = b(k, 1)
1994  tmp2 = b(k, 2)
1995  tmp3 = b(k, 3)
1996  tmp4 = b(k, 4)
1997  tmp5 = b(k, 5)
1998  tmp6 = b(k, 6)
1999  tmp7 = b(k, 7)
2000  do i=1,n1
2001  c(i, 1) = c(i, 1) + a(i,k) * tmp1
2002  c(i, 2) = c(i, 2) + a(i,k) * tmp2
2003  c(i, 3) = c(i, 3) + a(i,k) * tmp3
2004  c(i, 4) = c(i, 4) + a(i,k) * tmp4
2005  c(i, 5) = c(i, 5) + a(i,k) * tmp5
2006  c(i, 6) = c(i, 6) + a(i,k) * tmp6
2007  c(i, 7) = c(i, 7) + a(i,k) * tmp7
2008  enddo
2009 
2010  enddo
2011 
2012  return
2013  end subroutine mxmur3_7
2014  subroutine mxmur3_6(a,n1,b,n2,c,n3)
2015  use kinds, only : dp
2016  real(DP) :: a(n1,n2),b(n2,6),c(n1,6)
2017 
2018  do k=1,n2
2019  tmp1 = b(k, 1)
2020  tmp2 = b(k, 2)
2021  tmp3 = b(k, 3)
2022  tmp4 = b(k, 4)
2023  tmp5 = b(k, 5)
2024  tmp6 = b(k, 6)
2025  do i=1,n1
2026  c(i, 1) = c(i, 1) + a(i,k) * tmp1
2027  c(i, 2) = c(i, 2) + a(i,k) * tmp2
2028  c(i, 3) = c(i, 3) + a(i,k) * tmp3
2029  c(i, 4) = c(i, 4) + a(i,k) * tmp4
2030  c(i, 5) = c(i, 5) + a(i,k) * tmp5
2031  c(i, 6) = c(i, 6) + a(i,k) * tmp6
2032  enddo
2033 
2034  enddo
2035 
2036  return
2037  end subroutine mxmur3_6
2038  subroutine mxmur3_5(a,n1,b,n2,c,n3)
2039  use kinds, only : dp
2040  real(DP) :: a(n1,n2),b(n2,5),c(n1,5)
2041 
2042  do k=1,n2
2043  tmp1 = b(k, 1)
2044  tmp2 = b(k, 2)
2045  tmp3 = b(k, 3)
2046  tmp4 = b(k, 4)
2047  tmp5 = b(k, 5)
2048  do i=1,n1
2049  c(i, 1) = c(i, 1) + a(i,k) * tmp1
2050  c(i, 2) = c(i, 2) + a(i,k) * tmp2
2051  c(i, 3) = c(i, 3) + a(i,k) * tmp3
2052  c(i, 4) = c(i, 4) + a(i,k) * tmp4
2053  c(i, 5) = c(i, 5) + a(i,k) * tmp5
2054  enddo
2055 
2056  enddo
2057 
2058  return
2059  end subroutine mxmur3_5
2060  subroutine mxmur3_4(a,n1,b,n2,c,n3)
2061  use kinds, only : dp
2062  real(DP) :: a(n1,n2),b(n2,4),c(n1,4)
2063 
2064  do k=1,n2
2065  tmp1 = b(k, 1)
2066  tmp2 = b(k, 2)
2067  tmp3 = b(k, 3)
2068  tmp4 = b(k, 4)
2069  do i=1,n1
2070  c(i, 1) = c(i, 1) + a(i,k) * tmp1
2071  c(i, 2) = c(i, 2) + a(i,k) * tmp2
2072  c(i, 3) = c(i, 3) + a(i,k) * tmp3
2073  c(i, 4) = c(i, 4) + a(i,k) * tmp4
2074  enddo
2075 
2076  enddo
2077 
2078  return
2079  end subroutine mxmur3_4
2080  subroutine mxmur3_3(a,n1,b,n2,c,n3)
2081  use kinds, only : dp
2082  real(DP) :: a(n1,n2),b(n2,3),c(n1,3)
2083 
2084  do k=1,n2
2085  tmp1 = b(k, 1)
2086  tmp2 = b(k, 2)
2087  tmp3 = b(k, 3)
2088  do i=1,n1
2089  c(i, 1) = c(i, 1) + a(i,k) * tmp1
2090  c(i, 2) = c(i, 2) + a(i,k) * tmp2
2091  c(i, 3) = c(i, 3) + a(i,k) * tmp3
2092  enddo
2093 
2094  enddo
2095 
2096  return
2097  end subroutine mxmur3_3
2098  subroutine mxmur3_2(a,n1,b,n2,c,n3)
2099  use kinds, only : dp
2100  real(DP) :: a(n1,n2),b(n2,2),c(n1,2)
2101 
2102  do k=1,n2
2103  tmp1 = b(k, 1)
2104  tmp2 = b(k, 2)
2105  do i=1,n1
2106  c(i, 1) = c(i, 1) + a(i,k) * tmp1
2107  c(i, 2) = c(i, 2) + a(i,k) * tmp2
2108  enddo
2109 
2110  enddo
2111 
2112  return
2113  end subroutine mxmur3_2
2114  subroutine mxmur3_1(a,n1,b,n2,c,n3)
2115  use kinds, only : dp
2116  real(DP) :: a(n1,n2),b(n2,1),c(n1,1)
2117 
2118  do k=1,n2
2119  tmp1 = b(k, 1)
2120  do i=1,n1
2121  c(i, 1) = c(i, 1) + a(i,k) * tmp1
2122  enddo
2123  enddo
2124 
2125  return
2126  end subroutine mxmur3_1
2127 !----------------------------------------------------------------------
2128  subroutine mxmd(a,n1,b,n2,c,n3)
2129 
2130 ! Matrix-vector product routine.
2131 ! NOTE: Use assembly coded routine if available.
2132 
2133 !---------------------------------------------------------------------
2134  use kinds, only : dp
2135  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
2136  REAL(DP) :: ONE,ZERO,EPS
2137 
2138 
2139 
2140  one=1.0
2141  zero=0.0
2142  call dgemm( 'N','N',n1,n3,n2,one,a,n1,b,n2,zero,c,n1)
2143  return
2144  end subroutine mxmd
2145 !-----------------------------------------------------------------------
2146  subroutine mxmfb(a,n1,b,n2,c,n3)
2147 !-----------------------------------------------------------------------
2148 
2149 ! Matrix-vector product routine.
2150 ! NOTE: Use assembly coded routine if available.
2151 
2152 !----------------------------------------------------------------------
2153  use kinds, only : dp
2154  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
2155 
2156  integer :: wdsize
2157  save wdsize
2158  data wdsize/0/
2159 
2160 ! First call: determine word size for dgemm/sgemm discrimination, below.
2161 
2162  if (wdsize == 0) then
2163  one = 1.0
2164  eps = 1.e-12
2165  wdsize = 8
2166  if (one+eps == 1.0) wdsize = 4
2167  endif
2168 
2169  if (n2 <= 8) then
2170  if (n2 == 1) then
2171  call mxmfb_1(a,n1,b,n2,c,n3)
2172  elseif (n2 == 2) then
2173  call mxmfb_2(a,n1,b,n2,c,n3)
2174  elseif (n2 == 3) then
2175  call mxmfb_3(a,n1,b,n2,c,n3)
2176  elseif (n2 == 4) then
2177  call mxmfb_4(a,n1,b,n2,c,n3)
2178  elseif (n2 == 5) then
2179  call mxmfb_5(a,n1,b,n2,c,n3)
2180  elseif (n2 == 6) then
2181  call mxmfb_6(a,n1,b,n2,c,n3)
2182  elseif (n2 == 7) then
2183  call mxmfb_7(a,n1,b,n2,c,n3)
2184  else
2185  call mxmfb_8(a,n1,b,n2,c,n3)
2186  endif
2187  elseif (n2 <= 16) then
2188  if (n2 == 9) then
2189  call mxmfb_9(a,n1,b,n2,c,n3)
2190  elseif (n2 == 10) then
2191  call mxmfb_10(a,n1,b,n2,c,n3)
2192  elseif (n2 == 11) then
2193  call mxmfb_11(a,n1,b,n2,c,n3)
2194  elseif (n2 == 12) then
2195  call mxmfb_12(a,n1,b,n2,c,n3)
2196  elseif (n2 == 13) then
2197  call mxmfb_13(a,n1,b,n2,c,n3)
2198  elseif (n2 == 14) then
2199  call mxmfb_14(a,n1,b,n2,c,n3)
2200  elseif (n2 == 15) then
2201  call mxmfb_15(a,n1,b,n2,c,n3)
2202  else
2203  call mxmfb_16(a,n1,b,n2,c,n3)
2204  endif
2205  elseif (n2 <= 24) then
2206  if (n2 == 17) then
2207  call mxmfb_17(a,n1,b,n2,c,n3)
2208  elseif (n2 == 18) then
2209  call mxmfb_18(a,n1,b,n2,c,n3)
2210  elseif (n2 == 19) then
2211  call mxmfb_19(a,n1,b,n2,c,n3)
2212  elseif (n2 == 20) then
2213  call mxmfb_20(a,n1,b,n2,c,n3)
2214  elseif (n2 == 21) then
2215  call mxmfb_21(a,n1,b,n2,c,n3)
2216  elseif (n2 == 22) then
2217  call mxmfb_22(a,n1,b,n2,c,n3)
2218  elseif (n2 == 23) then
2219  call mxmfb_23(a,n1,b,n2,c,n3)
2220  elseif (n2 == 24) then
2221  call mxmfb_24(a,n1,b,n2,c,n3)
2222  endif
2223  else
2224 
2225  one=1.0
2226  zero=0.0
2227  if (wdsize == 4) then
2228  call sgemm( 'N','N',n1,n3,n2,one,a,n1,b,n2,zero,c,n1)
2229  else
2230  call dgemm( 'N','N',n1,n3,n2,one,a,n1,b,n2,zero,c,n1)
2231  endif
2232 
2233  endif
2234  return
2235  end subroutine mxmfb
2236 !-----------------------------------------------------------------------
2237  subroutine mxmfb_1(a,n1,b,n2,c,n3)
2238 
2239  use kinds, only : dp
2240  real(DP) :: a(n1,1),b(1,n3),c(n1,n3)
2241 
2242  do j=1,n3
2243  do i=1,n1
2244  c(i,j) = a(i,1)*b(1,j)
2245  enddo
2246  enddo
2247  return
2248  end subroutine mxmfb_1
2249 !-----------------------------------------------------------------------
2250  subroutine mxmfb_2(a,n1,b,n2,c,n3)
2251 
2252  use kinds, only : dp
2253  real(DP) :: a(n1,2),b(2,n3),c(n1,n3)
2254 
2255  do j=1,n3
2256  do i=1,n1
2257  c(i,j) = a(i,1)*b(1,j) &
2258  + a(i,2)*b(2,j)
2259  enddo
2260  enddo
2261  return
2262  end subroutine mxmfb_2
2263 !-----------------------------------------------------------------------
2264  subroutine mxmfb_3(a,n1,b,n2,c,n3)
2265 
2266  use kinds, only : dp
2267  real(DP) :: a(n1,3),b(3,n3),c(n1,n3)
2268 
2269  do j=1,n3
2270  do i=1,n1
2271  c(i,j) = a(i,1)*b(1,j) &
2272  + a(i,2)*b(2,j) &
2273  + a(i,3)*b(3,j)
2274  enddo
2275  enddo
2276  return
2277  end subroutine mxmfb_3
2278 !-----------------------------------------------------------------------
2279  subroutine mxmfb_4(a,n1,b,n2,c,n3)
2280 
2281  use kinds, only : dp
2282  real(DP) :: a(n1,4),b(4,n3),c(n1,n3)
2283 
2284  do j=1,n3
2285  do i=1,n1
2286  c(i,j) = a(i,1)*b(1,j) &
2287  + a(i,2)*b(2,j) &
2288  + a(i,3)*b(3,j) &
2289  + a(i,4)*b(4,j)
2290  enddo
2291  enddo
2292  return
2293  end subroutine mxmfb_4
2294 !-----------------------------------------------------------------------
2295  subroutine mxmfb_5(a,n1,b,n2,c,n3)
2296 
2297  use kinds, only : dp
2298  real(DP) :: a(n1,5),b(5,n3),c(n1,n3)
2299 
2300  do j=1,n3
2301  do i=1,n1
2302  c(i,j) = a(i,1)*b(1,j) &
2303  + a(i,2)*b(2,j) &
2304  + a(i,3)*b(3,j) &
2305  + a(i,4)*b(4,j) &
2306  + a(i,5)*b(5,j)
2307  enddo
2308  enddo
2309  return
2310  end subroutine mxmfb_5
2311 !-----------------------------------------------------------------------
2312  subroutine mxmfb_6(a,n1,b,n2,c,n3)
2313 
2314  use kinds, only : dp
2315  real(DP) :: a(n1,6),b(6,n3),c(n1,n3)
2316 
2317  do j=1,n3
2318  do i=1,n1
2319  c(i,j) = a(i,1)*b(1,j) &
2320  + a(i,2)*b(2,j) &
2321  + a(i,3)*b(3,j) &
2322  + a(i,4)*b(4,j) &
2323  + a(i,5)*b(5,j) &
2324  + a(i,6)*b(6,j)
2325  enddo
2326  enddo
2327  return
2328  end subroutine mxmfb_6
2329 !-----------------------------------------------------------------------
2330  subroutine mxmfb_7(a,n1,b,n2,c,n3)
2331 
2332  use kinds, only : dp
2333  real(DP) :: a(n1,7),b(7,n3),c(n1,n3)
2334 
2335  do j=1,n3
2336  do i=1,n1
2337  c(i,j) = a(i,1)*b(1,j) &
2338  + a(i,2)*b(2,j) &
2339  + a(i,3)*b(3,j) &
2340  + a(i,4)*b(4,j) &
2341  + a(i,5)*b(5,j) &
2342  + a(i,6)*b(6,j) &
2343  + a(i,7)*b(7,j)
2344  enddo
2345  enddo
2346  return
2347  end subroutine mxmfb_7
2348 !-----------------------------------------------------------------------
2349  subroutine mxmfb_8(a,n1,b,n2,c,n3)
2350 
2351  use kinds, only : dp
2352  real(DP) :: a(n1,8),b(8,n3),c(n1,n3)
2353 
2354  do j=1,n3
2355  do i=1,n1
2356  c(i,j) = a(i,1)*b(1,j) &
2357  + a(i,2)*b(2,j) &
2358  + a(i,3)*b(3,j) &
2359  + a(i,4)*b(4,j) &
2360  + a(i,5)*b(5,j) &
2361  + a(i,6)*b(6,j) &
2362  + a(i,7)*b(7,j) &
2363  + a(i,8)*b(8,j)
2364  enddo
2365  enddo
2366  return
2367  end subroutine mxmfb_8
2368 !-----------------------------------------------------------------------
2369  subroutine mxmfb_9(a,n1,b,n2,c,n3)
2370 
2371  use kinds, only : dp
2372  real(DP) :: a(n1,9),b(9,n3),c(n1,n3)
2373 
2374  do j=1,n3
2375  do i=1,n1
2376  c(i,j) = a(i,1)*b(1,j) &
2377  + a(i,2)*b(2,j) &
2378  + a(i,3)*b(3,j) &
2379  + a(i,4)*b(4,j) &
2380  + a(i,5)*b(5,j) &
2381  + a(i,6)*b(6,j) &
2382  + a(i,7)*b(7,j) &
2383  + a(i,8)*b(8,j) &
2384  + a(i,9)*b(9,j)
2385  enddo
2386  enddo
2387  return
2388  end subroutine mxmfb_9
2389 !-----------------------------------------------------------------------
2390  subroutine mxmfb_10(a,n1,b,n2,c,n3)
2391 
2392  use kinds, only : dp
2393  real(DP) :: a(n1,10),b(10,n3),c(n1,n3)
2394 
2395  do j=1,n3
2396  do i=1,n1
2397  c(i,j) = a(i,1)*b(1,j) &
2398  + a(i,2)*b(2,j) &
2399  + a(i,3)*b(3,j) &
2400  + a(i,4)*b(4,j) &
2401  + a(i,5)*b(5,j) &
2402  + a(i,6)*b(6,j) &
2403  + a(i,7)*b(7,j) &
2404  + a(i,8)*b(8,j) &
2405  + a(i,9)*b(9,j) &
2406  + a(i,10)*b(10,j)
2407  enddo
2408  enddo
2409  return
2410  end subroutine mxmfb_10
2411 !-----------------------------------------------------------------------
2412  subroutine mxmfb_11(a,n1,b,n2,c,n3)
2413 
2414  use kinds, only : dp
2415  real(DP) :: a(n1,11),b(11,n3),c(n1,n3)
2416 
2417  do j=1,n3
2418  do i=1,n1
2419  c(i,j) = a(i,1)*b(1,j) &
2420  + a(i,2)*b(2,j) &
2421  + a(i,3)*b(3,j) &
2422  + a(i,4)*b(4,j) &
2423  + a(i,5)*b(5,j) &
2424  + a(i,6)*b(6,j) &
2425  + a(i,7)*b(7,j) &
2426  + a(i,8)*b(8,j) &
2427  + a(i,9)*b(9,j) &
2428  + a(i,10)*b(10,j) &
2429  + a(i,11)*b(11,j)
2430  enddo
2431  enddo
2432  return
2433  end subroutine mxmfb_11
2434 !-----------------------------------------------------------------------
2435  subroutine mxmfb_12(a,n1,b,n2,c,n3)
2436 
2437  use kinds, only : dp
2438  real(DP) :: a(n1,12),b(12,n3),c(n1,n3)
2439 
2440  do j=1,n3
2441  do i=1,n1
2442  c(i,j) = a(i,1)*b(1,j) &
2443  + a(i,2)*b(2,j) &
2444  + a(i,3)*b(3,j) &
2445  + a(i,4)*b(4,j) &
2446  + a(i,5)*b(5,j) &
2447  + a(i,6)*b(6,j) &
2448  + a(i,7)*b(7,j) &
2449  + a(i,8)*b(8,j) &
2450  + a(i,9)*b(9,j) &
2451  + a(i,10)*b(10,j) &
2452  + a(i,11)*b(11,j) &
2453  + a(i,12)*b(12,j)
2454  enddo
2455  enddo
2456  return
2457  end subroutine mxmfb_12
2458 !-----------------------------------------------------------------------
2459  subroutine mxmfb_13(a,n1,b,n2,c,n3)
2460 
2461  use kinds, only : dp
2462  real(DP) :: a(n1,13),b(13,n3),c(n1,n3)
2463 
2464  do j=1,n3
2465  do i=1,n1
2466  c(i,j) = a(i,1)*b(1,j) &
2467  + a(i,2)*b(2,j) &
2468  + a(i,3)*b(3,j) &
2469  + a(i,4)*b(4,j) &
2470  + a(i,5)*b(5,j) &
2471  + a(i,6)*b(6,j) &
2472  + a(i,7)*b(7,j) &
2473  + a(i,8)*b(8,j) &
2474  + a(i,9)*b(9,j) &
2475  + a(i,10)*b(10,j) &
2476  + a(i,11)*b(11,j) &
2477  + a(i,12)*b(12,j) &
2478  + a(i,13)*b(13,j)
2479  enddo
2480  enddo
2481  return
2482  end subroutine mxmfb_13
2483 !-----------------------------------------------------------------------
2484  subroutine mxmfb_14(a,n1,b,n2,c,n3)
2485 
2486  use kinds, only : dp
2487  real(DP) :: a(n1,14),b(14,n3),c(n1,n3)
2488 
2489  do j=1,n3
2490  do i=1,n1
2491  c(i,j) = a(i,1)*b(1,j) &
2492  + a(i,2)*b(2,j) &
2493  + a(i,3)*b(3,j) &
2494  + a(i,4)*b(4,j) &
2495  + a(i,5)*b(5,j) &
2496  + a(i,6)*b(6,j) &
2497  + a(i,7)*b(7,j) &
2498  + a(i,8)*b(8,j) &
2499  + a(i,9)*b(9,j) &
2500  + a(i,10)*b(10,j) &
2501  + a(i,11)*b(11,j) &
2502  + a(i,12)*b(12,j) &
2503  + a(i,13)*b(13,j) &
2504  + a(i,14)*b(14,j)
2505  enddo
2506  enddo
2507  return
2508  end subroutine mxmfb_14
2509 !-----------------------------------------------------------------------
2510  subroutine mxmfb_15(a,n1,b,n2,c,n3)
2511 
2512  use kinds, only : dp
2513  real(DP) :: a(n1,15),b(15,n3),c(n1,n3)
2514 
2515  do j=1,n3
2516  do i=1,n1
2517  c(i,j) = a(i,1)*b(1,j) &
2518  + a(i,2)*b(2,j) &
2519  + a(i,3)*b(3,j) &
2520  + a(i,4)*b(4,j) &
2521  + a(i,5)*b(5,j) &
2522  + a(i,6)*b(6,j) &
2523  + a(i,7)*b(7,j) &
2524  + a(i,8)*b(8,j) &
2525  + a(i,9)*b(9,j) &
2526  + a(i,10)*b(10,j) &
2527  + a(i,11)*b(11,j) &
2528  + a(i,12)*b(12,j) &
2529  + a(i,13)*b(13,j) &
2530  + a(i,14)*b(14,j) &
2531  + a(i,15)*b(15,j)
2532  enddo
2533  enddo
2534  return
2535  end subroutine mxmfb_15
2536 !-----------------------------------------------------------------------
2537  subroutine mxmfb_16(a,n1,b,n2,c,n3)
2538 
2539  use kinds, only : dp
2540  real(DP) :: a(n1,16),b(16,n3),c(n1,n3)
2541 
2542  do j=1,n3
2543  do i=1,n1
2544  c(i,j) = a(i,1)*b(1,j) &
2545  + a(i,2)*b(2,j) &
2546  + a(i,3)*b(3,j) &
2547  + a(i,4)*b(4,j) &
2548  + a(i,5)*b(5,j) &
2549  + a(i,6)*b(6,j) &
2550  + a(i,7)*b(7,j) &
2551  + a(i,8)*b(8,j) &
2552  + a(i,9)*b(9,j) &
2553  + a(i,10)*b(10,j) &
2554  + a(i,11)*b(11,j) &
2555  + a(i,12)*b(12,j) &
2556  + a(i,13)*b(13,j) &
2557  + a(i,14)*b(14,j) &
2558  + a(i,15)*b(15,j) &
2559  + a(i,16)*b(16,j)
2560  enddo
2561  enddo
2562  return
2563  end subroutine mxmfb_16
2564 !-----------------------------------------------------------------------
2565  subroutine mxmfb_17(a,n1,b,n2,c,n3)
2566 
2567  use kinds, only : dp
2568  real(DP) :: a(n1,17),b(17,n3),c(n1,n3)
2569 
2570  do j=1,n3
2571  do i=1,n1
2572  c(i,j) = a(i,1)*b(1,j) &
2573  + a(i,2)*b(2,j) &
2574  + a(i,3)*b(3,j) &
2575  + a(i,4)*b(4,j) &
2576  + a(i,5)*b(5,j) &
2577  + a(i,6)*b(6,j) &
2578  + a(i,7)*b(7,j) &
2579  + a(i,8)*b(8,j) &
2580  + a(i,9)*b(9,j) &
2581  + a(i,10)*b(10,j) &
2582  + a(i,11)*b(11,j) &
2583  + a(i,12)*b(12,j) &
2584  + a(i,13)*b(13,j) &
2585  + a(i,14)*b(14,j) &
2586  + a(i,15)*b(15,j) &
2587  + a(i,16)*b(16,j) &
2588  + a(i,17)*b(17,j)
2589  enddo
2590  enddo
2591  return
2592  end subroutine mxmfb_17
2593 !-----------------------------------------------------------------------
2594  subroutine mxmfb_18(a,n1,b,n2,c,n3)
2595 
2596  use kinds, only : dp
2597  real(DP) :: a(n1,18),b(18,n3),c(n1,n3)
2598 
2599  do j=1,n3
2600  do i=1,n1
2601  c(i,j) = a(i,1)*b(1,j) &
2602  + a(i,2)*b(2,j) &
2603  + a(i,3)*b(3,j) &
2604  + a(i,4)*b(4,j) &
2605  + a(i,5)*b(5,j) &
2606  + a(i,6)*b(6,j) &
2607  + a(i,7)*b(7,j) &
2608  + a(i,8)*b(8,j) &
2609  + a(i,9)*b(9,j) &
2610  + a(i,10)*b(10,j) &
2611  + a(i,11)*b(11,j) &
2612  + a(i,12)*b(12,j) &
2613  + a(i,13)*b(13,j) &
2614  + a(i,14)*b(14,j) &
2615  + a(i,15)*b(15,j) &
2616  + a(i,16)*b(16,j) &
2617  + a(i,17)*b(17,j) &
2618  + a(i,18)*b(18,j)
2619  enddo
2620  enddo
2621  return
2622  end subroutine mxmfb_18
2623 !-----------------------------------------------------------------------
2624  subroutine mxmfb_19(a,n1,b,n2,c,n3)
2625 
2626  use kinds, only : dp
2627  real(DP) :: a(n1,19),b(19,n3),c(n1,n3)
2628 
2629  do j=1,n3
2630  do i=1,n1
2631  c(i,j) = a(i,1)*b(1,j) &
2632  + a(i,2)*b(2,j) &
2633  + a(i,3)*b(3,j) &
2634  + a(i,4)*b(4,j) &
2635  + a(i,5)*b(5,j) &
2636  + a(i,6)*b(6,j) &
2637  + a(i,7)*b(7,j) &
2638  + a(i,8)*b(8,j) &
2639  + a(i,9)*b(9,j) &
2640  + a(i,10)*b(10,j) &
2641  + a(i,11)*b(11,j) &
2642  + a(i,12)*b(12,j) &
2643  + a(i,13)*b(13,j) &
2644  + a(i,14)*b(14,j) &
2645  + a(i,15)*b(15,j) &
2646  + a(i,16)*b(16,j) &
2647  + a(i,17)*b(17,j) &
2648  + a(i,18)*b(18,j) &
2649  + a(i,19)*b(19,j)
2650  enddo
2651  enddo
2652  return
2653  end subroutine mxmfb_19
2654 !-----------------------------------------------------------------------
2655  subroutine mxmfb_20(a,n1,b,n2,c,n3)
2656 
2657  use kinds, only : dp
2658  real(DP) :: a(n1,20),b(20,n3),c(n1,n3)
2659 
2660  do j=1,n3
2661  do i=1,n1
2662  c(i,j) = a(i,1)*b(1,j) &
2663  + a(i,2)*b(2,j) &
2664  + a(i,3)*b(3,j) &
2665  + a(i,4)*b(4,j) &
2666  + a(i,5)*b(5,j) &
2667  + a(i,6)*b(6,j) &
2668  + a(i,7)*b(7,j) &
2669  + a(i,8)*b(8,j) &
2670  + a(i,9)*b(9,j) &
2671  + a(i,10)*b(10,j) &
2672  + a(i,11)*b(11,j) &
2673  + a(i,12)*b(12,j) &
2674  + a(i,13)*b(13,j) &
2675  + a(i,14)*b(14,j) &
2676  + a(i,15)*b(15,j) &
2677  + a(i,16)*b(16,j) &
2678  + a(i,17)*b(17,j) &
2679  + a(i,18)*b(18,j) &
2680  + a(i,19)*b(19,j) &
2681  + a(i,20)*b(20,j)
2682  enddo
2683  enddo
2684  return
2685  end subroutine mxmfb_20
2686 !-----------------------------------------------------------------------
2687  subroutine mxmfb_21(a,n1,b,n2,c,n3)
2688 
2689  use kinds, only : dp
2690  real(DP) :: a(n1,21),b(21,n3),c(n1,n3)
2691 
2692  do j=1,n3
2693  do i=1,n1
2694  c(i,j) = a(i,1)*b(1,j) &
2695  + a(i,2)*b(2,j) &
2696  + a(i,3)*b(3,j) &
2697  + a(i,4)*b(4,j) &
2698  + a(i,5)*b(5,j) &
2699  + a(i,6)*b(6,j) &
2700  + a(i,7)*b(7,j) &
2701  + a(i,8)*b(8,j) &
2702  + a(i,9)*b(9,j) &
2703  + a(i,10)*b(10,j) &
2704  + a(i,11)*b(11,j) &
2705  + a(i,12)*b(12,j) &
2706  + a(i,13)*b(13,j) &
2707  + a(i,14)*b(14,j) &
2708  + a(i,15)*b(15,j) &
2709  + a(i,16)*b(16,j) &
2710  + a(i,17)*b(17,j) &
2711  + a(i,18)*b(18,j) &
2712  + a(i,19)*b(19,j) &
2713  + a(i,20)*b(20,j) &
2714  + a(i,21)*b(21,j)
2715  enddo
2716  enddo
2717  return
2718  end subroutine mxmfb_21
2719 !-----------------------------------------------------------------------
2720  subroutine mxmfb_22(a,n1,b,n2,c,n3)
2721 
2722  use kinds, only : dp
2723  real(DP) :: a(n1,22),b(22,n3),c(n1,n3)
2724 
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  enddo
2750  enddo
2751  return
2752  end subroutine mxmfb_22
2753 !-----------------------------------------------------------------------
2754  subroutine mxmfb_23(a,n1,b,n2,c,n3)
2755 
2756  use kinds, only : dp
2757  real(DP) :: a(n1,23),b(23,n3),c(n1,n3)
2758 
2759  do j=1,n3
2760  do i=1,n1
2761  c(i,j) = a(i,1)*b(1,j) &
2762  + a(i,2)*b(2,j) &
2763  + a(i,3)*b(3,j) &
2764  + a(i,4)*b(4,j) &
2765  + a(i,5)*b(5,j) &
2766  + a(i,6)*b(6,j) &
2767  + a(i,7)*b(7,j) &
2768  + a(i,8)*b(8,j) &
2769  + a(i,9)*b(9,j) &
2770  + a(i,10)*b(10,j) &
2771  + a(i,11)*b(11,j) &
2772  + a(i,12)*b(12,j) &
2773  + a(i,13)*b(13,j) &
2774  + a(i,14)*b(14,j) &
2775  + a(i,15)*b(15,j) &
2776  + a(i,16)*b(16,j) &
2777  + a(i,17)*b(17,j) &
2778  + a(i,18)*b(18,j) &
2779  + a(i,19)*b(19,j) &
2780  + a(i,20)*b(20,j) &
2781  + a(i,21)*b(21,j) &
2782  + a(i,22)*b(22,j) &
2783  + a(i,23)*b(23,j)
2784  enddo
2785  enddo
2786  return
2787  end subroutine mxmfb_23
2788 !-----------------------------------------------------------------------
2789  subroutine mxmfb_24(a,n1,b,n2,c,n3)
2790 
2791  use kinds, only : dp
2792  real(DP) :: a(n1,24),b(24,n3),c(n1,n3)
2793 
2794  do j=1,n3
2795  do i=1,n1
2796  c(i,j) = a(i,1)*b(1,j) &
2797  + a(i,2)*b(2,j) &
2798  + a(i,3)*b(3,j) &
2799  + a(i,4)*b(4,j) &
2800  + a(i,5)*b(5,j) &
2801  + a(i,6)*b(6,j) &
2802  + a(i,7)*b(7,j) &
2803  + a(i,8)*b(8,j) &
2804  + a(i,9)*b(9,j) &
2805  + a(i,10)*b(10,j) &
2806  + a(i,11)*b(11,j) &
2807  + a(i,12)*b(12,j) &
2808  + a(i,13)*b(13,j) &
2809  + a(i,14)*b(14,j) &
2810  + a(i,15)*b(15,j) &
2811  + a(i,16)*b(16,j) &
2812  + a(i,17)*b(17,j) &
2813  + a(i,18)*b(18,j) &
2814  + a(i,19)*b(19,j) &
2815  + a(i,20)*b(20,j) &
2816  + a(i,21)*b(21,j) &
2817  + a(i,22)*b(22,j) &
2818  + a(i,23)*b(23,j) &
2819  + a(i,24)*b(24,j)
2820  enddo
2821  enddo
2822  return
2823  end subroutine mxmfb_24
2824 !-----------------------------------------------------------------------
2825  subroutine mxmf3(a,n1,b,n2,c,n3)
2826 !-----------------------------------------------------------------------
2827 
2828 ! Matrix-vector product routine.
2829 ! NOTE: Use assembly coded routine if available.
2830 
2831 !----------------------------------------------------------------------
2832  use kinds, only : dp
2833  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
2834 
2835  integer :: wdsize
2836  save wdsize
2837  data wdsize/0/
2838 
2839 ! First call: determine word size for dgemm/sgemm discrimination, below.
2840 
2841  if (wdsize == 0) then
2842  one = 1.0
2843  eps = 1.e-12
2844  wdsize = 8
2845  if (one+eps == 1.0) wdsize = 4
2846  endif
2847 
2848  if (n2 <= 8) then
2849  if (n2 == 1) then
2850  call mxmf3_1(a,n1,b,n2,c,n3)
2851  elseif (n2 == 2) then
2852  call mxmf3_2(a,n1,b,n2,c,n3)
2853  elseif (n2 == 3) then
2854  call mxmf3_3(a,n1,b,n2,c,n3)
2855  elseif (n2 == 4) then
2856  call mxmf3_4(a,n1,b,n2,c,n3)
2857  elseif (n2 == 5) then
2858  call mxmf3_5(a,n1,b,n2,c,n3)
2859  elseif (n2 == 6) then
2860  call mxmf3_6(a,n1,b,n2,c,n3)
2861  elseif (n2 == 7) then
2862  call mxmf3_7(a,n1,b,n2,c,n3)
2863  else
2864  call mxmf3_8(a,n1,b,n2,c,n3)
2865  endif
2866  elseif (n2 <= 16) then
2867  if (n2 == 9) then
2868  call mxmf3_9(a,n1,b,n2,c,n3)
2869  elseif (n2 == 10) then
2870  call mxmf3_10(a,n1,b,n2,c,n3)
2871  elseif (n2 == 11) then
2872  call mxmf3_11(a,n1,b,n2,c,n3)
2873  elseif (n2 == 12) then
2874  call mxmf3_12(a,n1,b,n2,c,n3)
2875  elseif (n2 == 13) then
2876  call mxmf3_13(a,n1,b,n2,c,n3)
2877  elseif (n2 == 14) then
2878  call mxmf3_14(a,n1,b,n2,c,n3)
2879  elseif (n2 == 15) then
2880  call mxmf3_15(a,n1,b,n2,c,n3)
2881  else
2882  call mxmf3_16(a,n1,b,n2,c,n3)
2883  endif
2884  elseif (n2 <= 24) then
2885  if (n2 == 17) then
2886  call mxmf3_17(a,n1,b,n2,c,n3)
2887  elseif (n2 == 18) then
2888  call mxmf3_18(a,n1,b,n2,c,n3)
2889  elseif (n2 == 19) then
2890  call mxmf3_19(a,n1,b,n2,c,n3)
2891  elseif (n2 == 20) then
2892  call mxmf3_20(a,n1,b,n2,c,n3)
2893  elseif (n2 == 21) then
2894  call mxmf3_21(a,n1,b,n2,c,n3)
2895  elseif (n2 == 22) then
2896  call mxmf3_22(a,n1,b,n2,c,n3)
2897  elseif (n2 == 23) then
2898  call mxmf3_23(a,n1,b,n2,c,n3)
2899  elseif (n2 == 24) then
2900  call mxmf3_24(a,n1,b,n2,c,n3)
2901  endif
2902  else
2903 
2904  one=1.0
2905  zero=0.0
2906  if (wdsize == 4) then
2907  call sgemm( 'N','N',n1,n3,n2,one,a,n1,b,n2,zero,c,n1)
2908  else
2909  call dgemm( 'N','N',n1,n3,n2,one,a,n1,b,n2,zero,c,n1)
2910  endif
2911 
2912  ! N0=N1*N3
2913  ! DO 10 I=1,N0
2914  ! C(I,1)=0.
2915  ! 10 CONTINUE
2916  ! DO 100 J=1,N3
2917  ! DO 100 K=1,N2
2918  ! BB=B(K,J)
2919  ! DO 100 I=1,N1
2920  ! C(I,J)=C(I,J)+A(I,K)*BB
2921  ! 100 CONTINUE
2922 
2923  endif
2924  return
2925  end subroutine mxmf3
2926 !-----------------------------------------------------------------------
2927  subroutine mxmf3_1(a,n1,b,n2,c,n3)
2928 
2929  use kinds, only : dp
2930  real(DP) :: a(n1,1),b(1,n3),c(n1,n3)
2931 
2932  do i=1,n1
2933  do j=1,n3
2934  c(i,j) = a(i,1)*b(1,j)
2935  enddo
2936  enddo
2937  return
2938  end subroutine mxmf3_1
2939 !-----------------------------------------------------------------------
2940  subroutine mxmf3_2(a,n1,b,n2,c,n3)
2941 
2942  use kinds, only : dp
2943  real(DP) :: a(n1,2),b(2,n3),c(n1,n3)
2944 
2945  do i=1,n1
2946  do j=1,n3
2947  c(i,j) = a(i,1)*b(1,j) &
2948  + a(i,2)*b(2,j)
2949  enddo
2950  enddo
2951  return
2952  end subroutine mxmf3_2
2953 !-----------------------------------------------------------------------
2954  subroutine mxmf3_3(a,n1,b,n2,c,n3)
2955 
2956  use kinds, only : dp
2957  real(DP) :: a(n1,3),b(3,n3),c(n1,n3)
2958 
2959  do i=1,n1
2960  do j=1,n3
2961  c(i,j) = a(i,1)*b(1,j) &
2962  + a(i,2)*b(2,j) &
2963  + a(i,3)*b(3,j)
2964  enddo
2965  enddo
2966  return
2967  end subroutine mxmf3_3
2968 !-----------------------------------------------------------------------
2969  subroutine mxmf3_4(a,n1,b,n2,c,n3)
2970 
2971  use kinds, only : dp
2972  real(DP) :: a(n1,4),b(4,n3),c(n1,n3)
2973 
2974  do i=1,n1
2975  do j=1,n3
2976  c(i,j) = a(i,1)*b(1,j) &
2977  + a(i,2)*b(2,j) &
2978  + a(i,3)*b(3,j) &
2979  + a(i,4)*b(4,j)
2980  enddo
2981  enddo
2982  return
2983  end subroutine mxmf3_4
2984 !-----------------------------------------------------------------------
2985  subroutine mxmf3_5(a,n1,b,n2,c,n3)
2986 
2987  use kinds, only : dp
2988  real(DP) :: a(n1,5),b(5,n3),c(n1,n3)
2989 
2990  do i=1,n1
2991  do j=1,n3
2992  c(i,j) = a(i,1)*b(1,j) &
2993  + a(i,2)*b(2,j) &
2994  + a(i,3)*b(3,j) &
2995  + a(i,4)*b(4,j) &
2996  + a(i,5)*b(5,j)
2997  enddo
2998  enddo
2999  return
3000  end subroutine mxmf3_5
3001 !-----------------------------------------------------------------------
3002  subroutine mxmf3_6(a,n1,b,n2,c,n3)
3003 
3004  use kinds, only : dp
3005  real(DP) :: a(n1,6),b(6,n3),c(n1,n3)
3006 
3007  do i=1,n1
3008  do j=1,n3
3009  c(i,j) = a(i,1)*b(1,j) &
3010  + a(i,2)*b(2,j) &
3011  + a(i,3)*b(3,j) &
3012  + a(i,4)*b(4,j) &
3013  + a(i,5)*b(5,j) &
3014  + a(i,6)*b(6,j)
3015  enddo
3016  enddo
3017  return
3018  end subroutine mxmf3_6
3019 !-----------------------------------------------------------------------
3020  subroutine mxmf3_7(a,n1,b,n2,c,n3)
3021 
3022  use kinds, only : dp
3023  real(DP) :: a(n1,7),b(7,n3),c(n1,n3)
3024 
3025  do i=1,n1
3026  do j=1,n3
3027  c(i,j) = a(i,1)*b(1,j) &
3028  + a(i,2)*b(2,j) &
3029  + a(i,3)*b(3,j) &
3030  + a(i,4)*b(4,j) &
3031  + a(i,5)*b(5,j) &
3032  + a(i,6)*b(6,j) &
3033  + a(i,7)*b(7,j)
3034  enddo
3035  enddo
3036  return
3037  end subroutine mxmf3_7
3038 !-----------------------------------------------------------------------
3039  subroutine mxmf3_8(a,n1,b,n2,c,n3)
3040 
3041  use kinds, only : dp
3042  real(DP) :: a(n1,8),b(8,n3),c(n1,n3)
3043 
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  enddo
3055  enddo
3056  return
3057  end subroutine mxmf3_8
3058 !-----------------------------------------------------------------------
3059  subroutine mxmf3_9(a,n1,b,n2,c,n3)
3060 
3061  use kinds, only : dp
3062  real(DP) :: a(n1,9),b(9,n3),c(n1,n3)
3063 
3064  do i=1,n1
3065  do j=1,n3
3066  c(i,j) = a(i,1)*b(1,j) &
3067  + a(i,2)*b(2,j) &
3068  + a(i,3)*b(3,j) &
3069  + a(i,4)*b(4,j) &
3070  + a(i,5)*b(5,j) &
3071  + a(i,6)*b(6,j) &
3072  + a(i,7)*b(7,j) &
3073  + a(i,8)*b(8,j) &
3074  + a(i,9)*b(9,j)
3075  enddo
3076  enddo
3077  return
3078  end subroutine mxmf3_9
3079 !-----------------------------------------------------------------------
3080  subroutine mxmf3_10(a,n1,b,n2,c,n3)
3081 
3082  use kinds, only : dp
3083  real(DP) :: a(n1,10),b(10,n3),c(n1,n3)
3084 
3085  do i=1,n1
3086  do j=1,n3
3087  c(i,j) = a(i,1)*b(1,j) &
3088  + a(i,2)*b(2,j) &
3089  + a(i,3)*b(3,j) &
3090  + a(i,4)*b(4,j) &
3091  + a(i,5)*b(5,j) &
3092  + a(i,6)*b(6,j) &
3093  + a(i,7)*b(7,j) &
3094  + a(i,8)*b(8,j) &
3095  + a(i,9)*b(9,j) &
3096  + a(i,10)*b(10,j)
3097  enddo
3098  enddo
3099  return
3100  end subroutine mxmf3_10
3101 !-----------------------------------------------------------------------
3102  subroutine mxmf3_11(a,n1,b,n2,c,n3)
3103 
3104  use kinds, only : dp
3105  real(DP) :: a(n1,11),b(11,n3),c(n1,n3)
3106 
3107  do i=1,n1
3108  do j=1,n3
3109  c(i,j) = a(i,1)*b(1,j) &
3110  + a(i,2)*b(2,j) &
3111  + a(i,3)*b(3,j) &
3112  + a(i,4)*b(4,j) &
3113  + a(i,5)*b(5,j) &
3114  + a(i,6)*b(6,j) &
3115  + a(i,7)*b(7,j) &
3116  + a(i,8)*b(8,j) &
3117  + a(i,9)*b(9,j) &
3118  + a(i,10)*b(10,j) &
3119  + a(i,11)*b(11,j)
3120  enddo
3121  enddo
3122  return
3123  end subroutine mxmf3_11
3124 !-----------------------------------------------------------------------
3125  subroutine mxmf3_12(a,n1,b,n2,c,n3)
3126 
3127  use kinds, only : dp
3128  real(DP) :: a(n1,12),b(12,n3),c(n1,n3)
3129 
3130  do i=1,n1
3131  do j=1,n3
3132  c(i,j) = a(i,1)*b(1,j) &
3133  + a(i,2)*b(2,j) &
3134  + a(i,3)*b(3,j) &
3135  + a(i,4)*b(4,j) &
3136  + a(i,5)*b(5,j) &
3137  + a(i,6)*b(6,j) &
3138  + a(i,7)*b(7,j) &
3139  + a(i,8)*b(8,j) &
3140  + a(i,9)*b(9,j) &
3141  + a(i,10)*b(10,j) &
3142  + a(i,11)*b(11,j) &
3143  + a(i,12)*b(12,j)
3144  enddo
3145  enddo
3146  return
3147  end subroutine mxmf3_12
3148 !-----------------------------------------------------------------------
3149  subroutine mxmf3_13(a,n1,b,n2,c,n3)
3150 
3151  use kinds, only : dp
3152  real(DP) :: a(n1,13),b(13,n3),c(n1,n3)
3153 
3154  do i=1,n1
3155  do j=1,n3
3156  c(i,j) = a(i,1)*b(1,j) &
3157  + a(i,2)*b(2,j) &
3158  + a(i,3)*b(3,j) &
3159  + a(i,4)*b(4,j) &
3160  + a(i,5)*b(5,j) &
3161  + a(i,6)*b(6,j) &
3162  + a(i,7)*b(7,j) &
3163  + a(i,8)*b(8,j) &
3164  + a(i,9)*b(9,j) &
3165  + a(i,10)*b(10,j) &
3166  + a(i,11)*b(11,j) &
3167  + a(i,12)*b(12,j) &
3168  + a(i,13)*b(13,j)
3169  enddo
3170  enddo
3171  return
3172  end subroutine mxmf3_13
3173 !-----------------------------------------------------------------------
3174  subroutine mxmf3_14(a,n1,b,n2,c,n3)
3175 
3176  use kinds, only : dp
3177  real(DP) :: a(n1,14),b(14,n3),c(n1,n3)
3178 
3179  do i=1,n1
3180  do j=1,n3
3181  c(i,j) = a(i,1)*b(1,j) &
3182  + a(i,2)*b(2,j) &
3183  + a(i,3)*b(3,j) &
3184  + a(i,4)*b(4,j) &
3185  + a(i,5)*b(5,j) &
3186  + a(i,6)*b(6,j) &
3187  + a(i,7)*b(7,j) &
3188  + a(i,8)*b(8,j) &
3189  + a(i,9)*b(9,j) &
3190  + a(i,10)*b(10,j) &
3191  + a(i,11)*b(11,j) &
3192  + a(i,12)*b(12,j) &
3193  + a(i,13)*b(13,j) &
3194  + a(i,14)*b(14,j)
3195  enddo
3196  enddo
3197  return
3198  end subroutine mxmf3_14
3199 !-----------------------------------------------------------------------
3200  subroutine mxmf3_15(a,n1,b,n2,c,n3)
3201 
3202  use kinds, only : dp
3203  real(DP) :: a(n1,15),b(15,n3),c(n1,n3)
3204 
3205  do i=1,n1
3206  do j=1,n3
3207  c(i,j) = a(i,1)*b(1,j) &
3208  + a(i,2)*b(2,j) &
3209  + a(i,3)*b(3,j) &
3210  + a(i,4)*b(4,j) &
3211  + a(i,5)*b(5,j) &
3212  + a(i,6)*b(6,j) &
3213  + a(i,7)*b(7,j) &
3214  + a(i,8)*b(8,j) &
3215  + a(i,9)*b(9,j) &
3216  + a(i,10)*b(10,j) &
3217  + a(i,11)*b(11,j) &
3218  + a(i,12)*b(12,j) &
3219  + a(i,13)*b(13,j) &
3220  + a(i,14)*b(14,j) &
3221  + a(i,15)*b(15,j)
3222  enddo
3223  enddo
3224  return
3225  end subroutine mxmf3_15
3226 !-----------------------------------------------------------------------
3227  subroutine mxmf3_16(a,n1,b,n2,c,n3)
3228 
3229  use kinds, only : dp
3230  real(DP) :: a(n1,16),b(16,n3),c(n1,n3)
3231 
3232  do i=1,n1
3233  do j=1,n3
3234  c(i,j) = a(i,1)*b(1,j) &
3235  + a(i,2)*b(2,j) &
3236  + a(i,3)*b(3,j) &
3237  + a(i,4)*b(4,j) &
3238  + a(i,5)*b(5,j) &
3239  + a(i,6)*b(6,j) &
3240  + a(i,7)*b(7,j) &
3241  + a(i,8)*b(8,j) &
3242  + a(i,9)*b(9,j) &
3243  + a(i,10)*b(10,j) &
3244  + a(i,11)*b(11,j) &
3245  + a(i,12)*b(12,j) &
3246  + a(i,13)*b(13,j) &
3247  + a(i,14)*b(14,j) &
3248  + a(i,15)*b(15,j) &
3249  + a(i,16)*b(16,j)
3250  enddo
3251  enddo
3252  return
3253  end subroutine mxmf3_16
3254 !-----------------------------------------------------------------------
3255  subroutine mxmf3_17(a,n1,b,n2,c,n3)
3256 
3257  use kinds, only : dp
3258  real(DP) :: a(n1,17),b(17,n3),c(n1,n3)
3259 
3260  do i=1,n1
3261  do j=1,n3
3262  c(i,j) = a(i,1)*b(1,j) &
3263  + a(i,2)*b(2,j) &
3264  + a(i,3)*b(3,j) &
3265  + a(i,4)*b(4,j) &
3266  + a(i,5)*b(5,j) &
3267  + a(i,6)*b(6,j) &
3268  + a(i,7)*b(7,j) &
3269  + a(i,8)*b(8,j) &
3270  + a(i,9)*b(9,j) &
3271  + a(i,10)*b(10,j) &
3272  + a(i,11)*b(11,j) &
3273  + a(i,12)*b(12,j) &
3274  + a(i,13)*b(13,j) &
3275  + a(i,14)*b(14,j) &
3276  + a(i,15)*b(15,j) &
3277  + a(i,16)*b(16,j) &
3278  + a(i,17)*b(17,j)
3279  enddo
3280  enddo
3281  return
3282  end subroutine mxmf3_17
3283 !-----------------------------------------------------------------------
3284  subroutine mxmf3_18(a,n1,b,n2,c,n3)
3285 
3286  use kinds, only : dp
3287  real(DP) :: a(n1,18),b(18,n3),c(n1,n3)
3288 
3289  do i=1,n1
3290  do j=1,n3
3291  c(i,j) = a(i,1)*b(1,j) &
3292  + a(i,2)*b(2,j) &
3293  + a(i,3)*b(3,j) &
3294  + a(i,4)*b(4,j) &
3295  + a(i,5)*b(5,j) &
3296  + a(i,6)*b(6,j) &
3297  + a(i,7)*b(7,j) &
3298  + a(i,8)*b(8,j) &
3299  + a(i,9)*b(9,j) &
3300  + a(i,10)*b(10,j) &
3301  + a(i,11)*b(11,j) &
3302  + a(i,12)*b(12,j) &
3303  + a(i,13)*b(13,j) &
3304  + a(i,14)*b(14,j) &
3305  + a(i,15)*b(15,j) &
3306  + a(i,16)*b(16,j) &
3307  + a(i,17)*b(17,j) &
3308  + a(i,18)*b(18,j)
3309  enddo
3310  enddo
3311  return
3312  end subroutine mxmf3_18
3313 !-----------------------------------------------------------------------
3314  subroutine mxmf3_19(a,n1,b,n2,c,n3)
3315 
3316  use kinds, only : dp
3317  real(DP) :: a(n1,19),b(19,n3),c(n1,n3)
3318 
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  enddo
3341  enddo
3342  return
3343  end subroutine mxmf3_19
3344 !-----------------------------------------------------------------------
3345  subroutine mxmf3_20(a,n1,b,n2,c,n3)
3346 
3347  use kinds, only : dp
3348  real(DP) :: a(n1,20),b(20,n3),c(n1,n3)
3349 
3350  do i=1,n1
3351  do j=1,n3
3352  c(i,j) = a(i,1)*b(1,j) &
3353  + a(i,2)*b(2,j) &
3354  + a(i,3)*b(3,j) &
3355  + a(i,4)*b(4,j) &
3356  + a(i,5)*b(5,j) &
3357  + a(i,6)*b(6,j) &
3358  + a(i,7)*b(7,j) &
3359  + a(i,8)*b(8,j) &
3360  + a(i,9)*b(9,j) &
3361  + a(i,10)*b(10,j) &
3362  + a(i,11)*b(11,j) &
3363  + a(i,12)*b(12,j) &
3364  + a(i,13)*b(13,j) &
3365  + a(i,14)*b(14,j) &
3366  + a(i,15)*b(15,j) &
3367  + a(i,16)*b(16,j) &
3368  + a(i,17)*b(17,j) &
3369  + a(i,18)*b(18,j) &
3370  + a(i,19)*b(19,j) &
3371  + a(i,20)*b(20,j)
3372  enddo
3373  enddo
3374  return
3375  end subroutine mxmf3_20
3376 !-----------------------------------------------------------------------
3377  subroutine mxmf3_21(a,n1,b,n2,c,n3)
3378 
3379  use kinds, only : dp
3380  real(DP) :: a(n1,21),b(21,n3),c(n1,n3)
3381 
3382  do i=1,n1
3383  do j=1,n3
3384  c(i,j) = a(i,1)*b(1,j) &
3385  + a(i,2)*b(2,j) &
3386  + a(i,3)*b(3,j) &
3387  + a(i,4)*b(4,j) &
3388  + a(i,5)*b(5,j) &
3389  + a(i,6)*b(6,j) &
3390  + a(i,7)*b(7,j) &
3391  + a(i,8)*b(8,j) &
3392  + a(i,9)*b(9,j) &
3393  + a(i,10)*b(10,j) &
3394  + a(i,11)*b(11,j) &
3395  + a(i,12)*b(12,j) &
3396  + a(i,13)*b(13,j) &
3397  + a(i,14)*b(14,j) &
3398  + a(i,15)*b(15,j) &
3399  + a(i,16)*b(16,j) &
3400  + a(i,17)*b(17,j) &
3401  + a(i,18)*b(18,j) &
3402  + a(i,19)*b(19,j) &
3403  + a(i,20)*b(20,j) &
3404  + a(i,21)*b(21,j)
3405  enddo
3406  enddo
3407  return
3408  end subroutine mxmf3_21
3409 !-----------------------------------------------------------------------
3410  subroutine mxmf3_22(a,n1,b,n2,c,n3)
3411 
3412  use kinds, only : dp
3413  real(DP) :: a(n1,22),b(22,n3),c(n1,n3)
3414 
3415  do i=1,n1
3416  do j=1,n3
3417  c(i,j) = a(i,1)*b(1,j) &
3418  + a(i,2)*b(2,j) &
3419  + a(i,3)*b(3,j) &
3420  + a(i,4)*b(4,j) &
3421  + a(i,5)*b(5,j) &
3422  + a(i,6)*b(6,j) &
3423  + a(i,7)*b(7,j) &
3424  + a(i,8)*b(8,j) &
3425  + a(i,9)*b(9,j) &
3426  + a(i,10)*b(10,j) &
3427  + a(i,11)*b(11,j) &
3428  + a(i,12)*b(12,j) &
3429  + a(i,13)*b(13,j) &
3430  + a(i,14)*b(14,j) &
3431  + a(i,15)*b(15,j) &
3432  + a(i,16)*b(16,j) &
3433  + a(i,17)*b(17,j) &
3434  + a(i,18)*b(18,j) &
3435  + a(i,19)*b(19,j) &
3436  + a(i,20)*b(20,j) &
3437  + a(i,21)*b(21,j) &
3438  + a(i,22)*b(22,j)
3439  enddo
3440  enddo
3441  return
3442  end subroutine mxmf3_22
3443 !-----------------------------------------------------------------------
3444  subroutine mxmf3_23(a,n1,b,n2,c,n3)
3445 
3446  use kinds, only : dp
3447  real(DP) :: a(n1,23),b(23,n3),c(n1,n3)
3448 
3449  do i=1,n1
3450  do j=1,n3
3451  c(i,j) = a(i,1)*b(1,j) &
3452  + a(i,2)*b(2,j) &
3453  + a(i,3)*b(3,j) &
3454  + a(i,4)*b(4,j) &
3455  + a(i,5)*b(5,j) &
3456  + a(i,6)*b(6,j) &
3457  + a(i,7)*b(7,j) &
3458  + a(i,8)*b(8,j) &
3459  + a(i,9)*b(9,j) &
3460  + a(i,10)*b(10,j) &
3461  + a(i,11)*b(11,j) &
3462  + a(i,12)*b(12,j) &
3463  + a(i,13)*b(13,j) &
3464  + a(i,14)*b(14,j) &
3465  + a(i,15)*b(15,j) &
3466  + a(i,16)*b(16,j) &
3467  + a(i,17)*b(17,j) &
3468  + a(i,18)*b(18,j) &
3469  + a(i,19)*b(19,j) &
3470  + a(i,20)*b(20,j) &
3471  + a(i,21)*b(21,j) &
3472  + a(i,22)*b(22,j) &
3473  + a(i,23)*b(23,j)
3474  enddo
3475  enddo
3476  return
3477  end subroutine mxmf3_23
3478 !-----------------------------------------------------------------------
3479  subroutine mxmf3_24(a,n1,b,n2,c,n3)
3480 
3481  use kinds, only : dp
3482  real(DP) :: a(n1,24),b(24,n3),c(n1,n3)
3483 
3484  do i=1,n1
3485  do j=1,n3
3486  c(i,j) = a(i,1)*b(1,j) &
3487  + a(i,2)*b(2,j) &
3488  + a(i,3)*b(3,j) &
3489  + a(i,4)*b(4,j) &
3490  + a(i,5)*b(5,j) &
3491  + a(i,6)*b(6,j) &
3492  + a(i,7)*b(7,j) &
3493  + a(i,8)*b(8,j) &
3494  + a(i,9)*b(9,j) &
3495  + a(i,10)*b(10,j) &
3496  + a(i,11)*b(11,j) &
3497  + a(i,12)*b(12,j) &
3498  + a(i,13)*b(13,j) &
3499  + a(i,14)*b(14,j) &
3500  + a(i,15)*b(15,j) &
3501  + a(i,16)*b(16,j) &
3502  + a(i,17)*b(17,j) &
3503  + a(i,18)*b(18,j) &
3504  + a(i,19)*b(19,j) &
3505  + a(i,20)*b(20,j) &
3506  + a(i,21)*b(21,j) &
3507  + a(i,22)*b(22,j) &
3508  + a(i,23)*b(23,j) &
3509  + a(i,24)*b(24,j)
3510  enddo
3511  enddo
3512  return
3513  end subroutine mxmf3_24
3514 !-----------------------------------------------------------------------
3515  subroutine mxm44(a,n1,b,n2,c,n3)
3516 !-----------------------------------------------------------------------
3517 
3518 ! NOTE -- this code has been set up with the "mxmf3" routine
3519 ! referenced in memtime.f. On most machines, the f2
3520 ! and f3 versions give the same performance (f2 is the
3521 ! nekton standard). On the t3e, f3 is noticeably faster.
3522 ! pff 10/5/98
3523 
3524 
3525 ! Matrix-vector product routine.
3526 ! NOTE: Use assembly coded routine if available.
3527 
3528 !----------------------------------------------------------------------
3529  use kinds, only : dp
3530  REAL(DP) :: A(n1,n2),B(n2,n3),C(n1,n3)
3531 
3532  if (n2 == 1) then
3533  call mxm44_2_t(a,n1,b,2,c,n3)
3534  elseif (n2 == 2) then
3535  call mxm44_2_t(a,n1,b,n2,c,n3)
3536  else
3537  call mxm44_0_t(a,n1,b,n2,c,n3)
3538  endif
3539 
3540  return
3541  end subroutine mxm44
3542 
3543 !-----------------------------------------------------------------------
3544  subroutine mxm44_0_t(a, m, b, k, c, n)
3545 ! subroutine matmul44(m, n, k, a, lda, b, ldb, c, ldc)
3546  use kinds, only : dp
3547 ! real*8 a(lda,k), b(ldb,n), c(ldc,n)
3548  real(DP) :: a(m,k), b(k,n), c(m,n)
3549  real(DP) :: s11, s12, s13, s14, s21, s22, s23, s24
3550  real(DP) :: s31, s32, s33, s34, s41, s42, s43, s44
3551 
3552 ! matrix multiply with a 4x4 pencil
3553 
3554 
3555  mresid = iand(m,3)
3556  nresid = iand(n,3)
3557  m1 = m - mresid + 1
3558  n1 = n - nresid + 1
3559 
3560  do i=1,m-mresid,4
3561  do j=1,n-nresid,4
3562  s11 = 0.0d0
3563  s21 = 0.0d0
3564  s31 = 0.0d0
3565  s41 = 0.0d0
3566  s12 = 0.0d0
3567  s22 = 0.0d0
3568  s32 = 0.0d0
3569  s42 = 0.0d0
3570  s13 = 0.0d0
3571  s23 = 0.0d0
3572  s33 = 0.0d0
3573  s43 = 0.0d0
3574  s14 = 0.0d0
3575  s24 = 0.0d0
3576  s34 = 0.0d0
3577  s44 = 0.0d0
3578  do l=1,k
3579  s11 = s11 + a(i,l)*b(l,j)
3580  s12 = s12 + a(i,l)*b(l,j+1)
3581  s13 = s13 + a(i,l)*b(l,j+2)
3582  s14 = s14 + a(i,l)*b(l,j+3)
3583 
3584  s21 = s21 + a(i+1,l)*b(l,j)
3585  s22 = s22 + a(i+1,l)*b(l,j+1)
3586  s23 = s23 + a(i+1,l)*b(l,j+2)
3587  s24 = s24 + a(i+1,l)*b(l,j+3)
3588 
3589  s31 = s31 + a(i+2,l)*b(l,j)
3590  s32 = s32 + a(i+2,l)*b(l,j+1)
3591  s33 = s33 + a(i+2,l)*b(l,j+2)
3592  s34 = s34 + a(i+2,l)*b(l,j+3)
3593 
3594  s41 = s41 + a(i+3,l)*b(l,j)
3595  s42 = s42 + a(i+3,l)*b(l,j+1)
3596  s43 = s43 + a(i+3,l)*b(l,j+2)
3597  s44 = s44 + a(i+3,l)*b(l,j+3)
3598  enddo
3599  c(i,j) = s11
3600  c(i,j+1) = s12
3601  c(i,j+2) = s13
3602  c(i,j+3) = s14
3603 
3604  c(i+1,j) = s21
3605  c(i+2,j) = s31
3606  c(i+3,j) = s41
3607 
3608  c(i+1,j+1) = s22
3609  c(i+2,j+1) = s32
3610  c(i+3,j+1) = s42
3611 
3612  c(i+1,j+2) = s23
3613  c(i+2,j+2) = s33
3614  c(i+3,j+2) = s43
3615 
3616  c(i+1,j+3) = s24
3617  c(i+2,j+3) = s34
3618  c(i+3,j+3) = s44
3619  enddo
3620  ! Residual when n is not multiple of 4
3621  if (nresid /= 0) then
3622  if (nresid == 1) then
3623  s11 = 0.0d0
3624  s21 = 0.0d0
3625  s31 = 0.0d0
3626  s41 = 0.0d0
3627  do l=1,k
3628  s11 = s11 + a(i,l)*b(l,n)
3629  s21 = s21 + a(i+1,l)*b(l,n)
3630  s31 = s31 + a(i+2,l)*b(l,n)
3631  s41 = s41 + a(i+3,l)*b(l,n)
3632  enddo
3633  c(i,n) = s11
3634  c(i+1,n) = s21
3635  c(i+2,n) = s31
3636  c(i+3,n) = s41
3637  elseif (nresid == 2) then
3638  s11 = 0.0d0
3639  s21 = 0.0d0
3640  s31 = 0.0d0
3641  s41 = 0.0d0
3642  s12 = 0.0d0
3643  s22 = 0.0d0
3644  s32 = 0.0d0
3645  s42 = 0.0d0
3646  do l=1,k
3647  s11 = s11 + a(i,l)*b(l,j)
3648  s12 = s12 + a(i,l)*b(l,j+1)
3649 
3650  s21 = s21 + a(i+1,l)*b(l,j)
3651  s22 = s22 + a(i+1,l)*b(l,j+1)
3652 
3653  s31 = s31 + a(i+2,l)*b(l,j)
3654  s32 = s32 + a(i+2,l)*b(l,j+1)
3655 
3656  s41 = s41 + a(i+3,l)*b(l,j)
3657  s42 = s42 + a(i+3,l)*b(l,j+1)
3658  enddo
3659  c(i,j) = s11
3660  c(i,j+1) = s12
3661 
3662  c(i+1,j) = s21
3663  c(i+2,j) = s31
3664  c(i+3,j) = s41
3665 
3666  c(i+1,j+1) = s22
3667  c(i+2,j+1) = s32
3668  c(i+3,j+1) = s42
3669  else
3670  s11 = 0.0d0
3671  s21 = 0.0d0
3672  s31 = 0.0d0
3673  s41 = 0.0d0
3674  s12 = 0.0d0
3675  s22 = 0.0d0
3676  s32 = 0.0d0
3677  s42 = 0.0d0
3678  s13 = 0.0d0
3679  s23 = 0.0d0
3680  s33 = 0.0d0
3681  s43 = 0.0d0
3682  do l=1,k
3683  s11 = s11 + a(i,l)*b(l,j)
3684  s12 = s12 + a(i,l)*b(l,j+1)
3685  s13 = s13 + a(i,l)*b(l,j+2)
3686 
3687  s21 = s21 + a(i+1,l)*b(l,j)
3688  s22 = s22 + a(i+1,l)*b(l,j+1)
3689  s23 = s23 + a(i+1,l)*b(l,j+2)
3690 
3691  s31 = s31 + a(i+2,l)*b(l,j)
3692  s32 = s32 + a(i+2,l)*b(l,j+1)
3693  s33 = s33 + a(i+2,l)*b(l,j+2)
3694 
3695  s41 = s41 + a(i+3,l)*b(l,j)
3696  s42 = s42 + a(i+3,l)*b(l,j+1)
3697  s43 = s43 + a(i+3,l)*b(l,j+2)
3698  enddo
3699  c(i,j) = s11
3700  c(i+1,j) = s21
3701  c(i+2,j) = s31
3702  c(i+3,j) = s41
3703  c(i,j+1) = s12
3704  c(i+1,j+1) = s22
3705  c(i+2,j+1) = s32
3706  c(i+3,j+1) = s42
3707  c(i,j+2) = s13
3708  c(i+1,j+2) = s23
3709  c(i+2,j+2) = s33
3710  c(i+3,j+2) = s43
3711  endif
3712  endif
3713  enddo
3714 
3715 ! Residual when m is not multiple of 4
3716  if (mresid == 0) then
3717  return
3718  elseif (mresid == 1) then
3719  do j=1,n-nresid,4
3720  s11 = 0.0d0
3721  s12 = 0.0d0
3722  s13 = 0.0d0
3723  s14 = 0.0d0
3724  do l=1,k
3725  s11 = s11 + a(m,l)*b(l,j)
3726  s12 = s12 + a(m,l)*b(l,j+1)
3727  s13 = s13 + a(m,l)*b(l,j+2)
3728  s14 = s14 + a(m,l)*b(l,j+3)
3729  enddo
3730  c(m,j) = s11
3731  c(m,j+1) = s12
3732  c(m,j+2) = s13
3733  c(m,j+3) = s14
3734  enddo
3735  ! mresid is 1, check nresid
3736  if (nresid == 0) then
3737  return
3738  elseif (nresid == 1) then
3739  s11 = 0.0d0
3740  do l=1,k
3741  s11 = s11 + a(m,l)*b(l,n)
3742  enddo
3743  c(m,n) = s11
3744  return
3745  elseif (nresid == 2) then
3746  s11 = 0.0d0
3747  s12 = 0.0d0
3748  do l=1,k
3749  s11 = s11 + a(m,l)*b(l,n-1)
3750  s12 = s12 + a(m,l)*b(l,n)
3751  enddo
3752  c(m,n-1) = s11
3753  c(m,n) = s12
3754  return
3755  else
3756  s11 = 0.0d0
3757  s12 = 0.0d0
3758  s13 = 0.0d0
3759  do l=1,k
3760  s11 = s11 + a(m,l)*b(l,n-2)
3761  s12 = s12 + a(m,l)*b(l,n-1)
3762  s13 = s13 + a(m,l)*b(l,n)
3763  enddo
3764  c(m,n-2) = s11
3765  c(m,n-1) = s12
3766  c(m,n) = s13
3767  return
3768  endif
3769  elseif (mresid == 2) then
3770  do j=1,n-nresid,4
3771  s11 = 0.0d0
3772  s12 = 0.0d0
3773  s13 = 0.0d0
3774  s14 = 0.0d0
3775  s21 = 0.0d0
3776  s22 = 0.0d0
3777  s23 = 0.0d0
3778  s24 = 0.0d0
3779  do l=1,k
3780  s11 = s11 + a(m-1,l)*b(l,j)
3781  s12 = s12 + a(m-1,l)*b(l,j+1)
3782  s13 = s13 + a(m-1,l)*b(l,j+2)
3783  s14 = s14 + a(m-1,l)*b(l,j+3)
3784 
3785  s21 = s21 + a(m,l)*b(l,j)
3786  s22 = s22 + a(m,l)*b(l,j+1)
3787  s23 = s23 + a(m,l)*b(l,j+2)
3788  s24 = s24 + a(m,l)*b(l,j+3)
3789  enddo
3790  c(m-1,j) = s11
3791  c(m-1,j+1) = s12
3792  c(m-1,j+2) = s13
3793  c(m-1,j+3) = s14
3794  c(m,j) = s21
3795  c(m,j+1) = s22
3796  c(m,j+2) = s23
3797  c(m,j+3) = s24
3798  enddo
3799  ! mresid is 2, check nresid
3800  if (nresid == 0) then
3801  return
3802  elseif (nresid == 1) then
3803  s11 = 0.0d0
3804  s21 = 0.0d0
3805  do l=1,k
3806  s11 = s11 + a(m-1,l)*b(l,n)
3807  s21 = s21 + a(m,l)*b(l,n)
3808  enddo
3809  c(m-1,n) = s11
3810  c(m,n) = s21
3811  return
3812  elseif (nresid == 2) then
3813  s11 = 0.0d0
3814  s21 = 0.0d0
3815  s12 = 0.0d0
3816  s22 = 0.0d0
3817  do l=1,k
3818  s11 = s11 + a(m-1,l)*b(l,n-1)
3819  s12 = s12 + a(m-1,l)*b(l,n)
3820  s21 = s21 + a(m,l)*b(l,n-1)
3821  s22 = s22 + a(m,l)*b(l,n)
3822  enddo
3823  c(m-1,n-1) = s11
3824  c(m-1,n) = s12
3825  c(m,n-1) = s21
3826  c(m,n) = s22
3827  return
3828  else
3829  s11 = 0.0d0
3830  s21 = 0.0d0
3831  s12 = 0.0d0
3832  s22 = 0.0d0
3833  s13 = 0.0d0
3834  s23 = 0.0d0
3835  do l=1,k
3836  s11 = s11 + a(m-1,l)*b(l,n-2)
3837  s12 = s12 + a(m-1,l)*b(l,n-1)
3838  s13 = s13 + a(m-1,l)*b(l,n)
3839  s21 = s21 + a(m,l)*b(l,n-2)
3840  s22 = s22 + a(m,l)*b(l,n-1)
3841  s23 = s23 + a(m,l)*b(l,n)
3842  enddo
3843  c(m-1,n-2) = s11
3844  c(m-1,n-1) = s12
3845  c(m-1,n) = s13
3846  c(m,n-2) = s21
3847  c(m,n-1) = s22
3848  c(m,n) = s23
3849  return
3850  endif
3851  else
3852  ! mresid is 3
3853  do j=1,n-nresid,4
3854  s11 = 0.0d0
3855  s21 = 0.0d0
3856  s31 = 0.0d0
3857 
3858  s12 = 0.0d0
3859  s22 = 0.0d0
3860  s32 = 0.0d0
3861 
3862  s13 = 0.0d0
3863  s23 = 0.0d0
3864  s33 = 0.0d0
3865 
3866  s14 = 0.0d0
3867  s24 = 0.0d0
3868  s34 = 0.0d0
3869 
3870  do l=1,k
3871  s11 = s11 + a(m-2,l)*b(l,j)
3872  s12 = s12 + a(m-2,l)*b(l,j+1)
3873  s13 = s13 + a(m-2,l)*b(l,j+2)
3874  s14 = s14 + a(m-2,l)*b(l,j+3)
3875 
3876  s21 = s21 + a(m-1,l)*b(l,j)
3877  s22 = s22 + a(m-1,l)*b(l,j+1)
3878  s23 = s23 + a(m-1,l)*b(l,j+2)
3879  s24 = s24 + a(m-1,l)*b(l,j+3)
3880 
3881  s31 = s31 + a(m,l)*b(l,j)
3882  s32 = s32 + a(m,l)*b(l,j+1)
3883  s33 = s33 + a(m,l)*b(l,j+2)
3884  s34 = s34 + a(m,l)*b(l,j+3)
3885  enddo
3886  c(m-2,j) = s11
3887  c(m-2,j+1) = s12
3888  c(m-2,j+2) = s13
3889  c(m-2,j+3) = s14
3890 
3891  c(m-1,j) = s21
3892  c(m-1,j+1) = s22
3893  c(m-1,j+2) = s23
3894  c(m-1,j+3) = s24
3895 
3896  c(m,j) = s31
3897  c(m,j+1) = s32
3898  c(m,j+2) = s33
3899  c(m,j+3) = s34
3900  enddo
3901  ! mresid is 3, check nresid
3902  if (nresid == 0) then
3903  return
3904  elseif (nresid == 1) then
3905  s11 = 0.0d0
3906  s21 = 0.0d0
3907  s31 = 0.0d0
3908  do l=1,k
3909  s11 = s11 + a(m-2,l)*b(l,n)
3910  s21 = s21 + a(m-1,l)*b(l,n)
3911  s31 = s31 + a(m,l)*b(l,n)
3912  enddo
3913  c(m-2,n) = s11
3914  c(m-1,n) = s21
3915  c(m,n) = s31
3916  return
3917  elseif (nresid == 2) then
3918  s11 = 0.0d0
3919  s21 = 0.0d0
3920  s31 = 0.0d0
3921  s12 = 0.0d0
3922  s22 = 0.0d0
3923  s32 = 0.0d0
3924  do l=1,k
3925  s11 = s11 + a(m-2,l)*b(l,n-1)
3926  s12 = s12 + a(m-2,l)*b(l,n)
3927  s21 = s21 + a(m-1,l)*b(l,n-1)
3928  s22 = s22 + a(m-1,l)*b(l,n)
3929  s31 = s31 + a(m,l)*b(l,n-1)
3930  s32 = s32 + a(m,l)*b(l,n)
3931  enddo
3932  c(m-2,n-1) = s11
3933  c(m-2,n) = s12
3934  c(m-1,n-1) = s21
3935  c(m-1,n) = s22
3936  c(m,n-1) = s31
3937  c(m,n) = s32
3938  return
3939  else
3940  s11 = 0.0d0
3941  s21 = 0.0d0
3942  s31 = 0.0d0
3943  s12 = 0.0d0
3944  s22 = 0.0d0
3945  s32 = 0.0d0
3946  s13 = 0.0d0
3947  s23 = 0.0d0
3948  s33 = 0.0d0
3949  do l=1,k
3950  s11 = s11 + a(m-2,l)*b(l,n-2)
3951  s12 = s12 + a(m-2,l)*b(l,n-1)
3952  s13 = s13 + a(m-2,l)*b(l,n)
3953  s21 = s21 + a(m-1,l)*b(l,n-2)
3954  s22 = s22 + a(m-1,l)*b(l,n-1)
3955  s23 = s23 + a(m-1,l)*b(l,n)
3956  s31 = s31 + a(m,l)*b(l,n-2)
3957  s32 = s32 + a(m,l)*b(l,n-1)
3958  s33 = s33 + a(m,l)*b(l,n)
3959  enddo
3960  c(m-2,n-2) = s11
3961  c(m-2,n-1) = s12
3962  c(m-2,n) = s13
3963  c(m-1,n-2) = s21
3964  c(m-1,n-1) = s22
3965  c(m-1,n) = s23
3966  c(m,n-2) = s31
3967  c(m,n-1) = s32
3968  c(m,n) = s33
3969  return
3970  endif
3971  endif
3972 
3973  return
3974  end subroutine mxm44_0_t
3975 !-----------------------------------------------------------------------
3976  subroutine mxm44_2_t(a, m, b, k, c, n)
3977  use kinds, only : dp
3978  real(DP) :: a(m,2), b(2,n), c(m,n)
3979 
3980  nresid = iand(n,3)
3981  n1 = n - nresid + 1
3982 
3983  do j=1,n-nresid,4
3984  do i=1,m
3985  c(i,j) = a(i,1)*b(1,j) &
3986  + a(i,2)*b(2,j)
3987  c(i,j+1) = a(i,1)*b(1,j+1) &
3988  + a(i,2)*b(2,j+1)
3989  c(i,j+2) = a(i,1)*b(1,j+2) &
3990  + a(i,2)*b(2,j+2)
3991  c(i,j+3) = a(i,1)*b(1,j+3) &
3992  + a(i,2)*b(2,j+3)
3993  enddo
3994  enddo
3995  if (nresid == 0) then
3996  return
3997  elseif (nresid == 1) then
3998  do i=1,m
3999  c(i,n) = a(i,1)*b(1,n) &
4000  + a(i,2)*b(2,n)
4001  enddo
4002  elseif (nresid == 2) then
4003  do i=1,m
4004  c(i,n-1) = a(i,1)*b(1,n-1) &
4005  + a(i,2)*b(2,n-1)
4006  c(i,n) = a(i,1)*b(1,n) &
4007  + a(i,2)*b(2,n)
4008  enddo
4009  else
4010  do i=1,m
4011  c(i,n-2) = a(i,1)*b(1,n-2) &
4012  + a(i,2)*b(2,n-2)
4013  c(i,n-1) = a(i,1)*b(1,n-1) &
4014  + a(i,2)*b(2,n-1)
4015  c(i,n) = a(i,1)*b(1,n) &
4016  + a(i,2)*b(2,n)
4017  enddo
4018  endif
4019 
4020  return
4021  end subroutine mxm44_2_t
4022 !-----------------------------------------------------------------------
mxmur2_11
subroutine mxmur2_11(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1455
mxmf3_11
subroutine mxmf3_11(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3102
mxmur2_1
subroutine mxmur2_1(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1290
mxmf2
subroutine mxmf2(A, N1, B, N2, C, N3)
unrolled loop version
Definition: mxm_std.F90:2
mxf7
subroutine mxf7(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:166
mxmfb_12
subroutine mxmfb_12(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2435
mxmur2_3
subroutine mxmur2_3(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1315
mxmfb_6
subroutine mxmfb_6(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2312
n
n
Definition: xxt_test.m:73
mxmur2_9
subroutine mxmur2_9(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1414
mxmf3_23
subroutine mxmf3_23(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3444
mxf14
subroutine mxf14(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:320
mxmf3_24
subroutine mxmf3_24(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3479
mxf9
subroutine mxf9(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:205
madd
subroutine madd(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1215
initab
subroutine initab(a, b, n)
Definition: mxm_std.F90:1145
mxmur3_2
subroutine mxmur3_2(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2098
mxmur3_16
subroutine mxmur3_16(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1664
mxmf3_6
subroutine mxmf3_6(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3002
mxmur2_16
subroutine mxmur2_16(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1575
mxmur2_13
subroutine mxmur2_13(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1500
mxmfb_16
subroutine mxmfb_16(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2537
mxm44_0
subroutine mxm44_0(a, m, b, k, c, n)
Definition: mxm_std.F90:668
mxmfb_11
subroutine mxmfb_11(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2412
mxmf3
subroutine mxmf3(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2825
mxmur2_8
subroutine mxmur2_8(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1395
mxmur3_4
subroutine mxmur3_4(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2060
mxmf3_15
subroutine mxmf3_15(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3200
mxf23
subroutine mxf23(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:597
mxmfb_7
subroutine mxmfb_7(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2330
mxmfb_13
subroutine mxmfb_13(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2459
mxmfb
subroutine mxmfb(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2146
mxf11
subroutine mxf11(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:248
mxm44_0_t
subroutine mxm44_0_t(a, m, b, k, c, n)
Definition: mxm_std.F90:3544
mxms
subroutine mxms(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1161
mxmur3_6
subroutine mxmur3_6(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2014
mxmfb_24
subroutine mxmfb_24(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2789
mxmur3_14
subroutine mxmur3_14(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1750
mxmfb_1
subroutine mxmfb_1(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2237
mxmf3_7
subroutine mxmf3_7(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3020
mxmf3_20
subroutine mxmf3_20(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3345
mxf1
subroutine mxf1(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:73
mxmur3_12
subroutine mxmur3_12(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1828
mxmf3_12
subroutine mxmf3_12(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3125
mxmd
subroutine mxmd(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2128
mxmfb_18
subroutine mxmfb_18(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2594
mxmfb_23
subroutine mxmfb_23(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2754
mxf19
subroutine mxf19(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:467
mxmfb_4
subroutine mxmfb_4(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2279
mxmur2_2
subroutine mxmur2_2(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1302
mxf21
subroutine mxf21(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:530
mxmfb_8
subroutine mxmfb_8(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2349
mxmfb_22
subroutine mxmfb_22(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2720
mxf6
subroutine mxf6(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:148
mxmf3_4
subroutine mxmf3_4(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2969
mxmfb_14
subroutine mxmfb_14(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2484
mxf20
subroutine mxf20(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:498
mxmf3_19
subroutine mxmf3_19(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3314
mxmfb_20
subroutine mxmfb_20(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2655
mxf10
subroutine mxf10(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:226
mxmur2
subroutine mxmur2(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1229
mxmf3_16
subroutine mxmf3_16(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3227
mxf15
subroutine mxf15(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:346
mxmur2_6
subroutine mxmur2_6(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1360
mxmf3_5
subroutine mxmf3_5(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2985
mxmur2_4
subroutine mxmur2_4(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1329
mxf13
subroutine mxf13(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:295
mxmfb_3
subroutine mxmfb_3(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2264
mxmf3_21
subroutine mxmf3_21(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3377
mxmur3
subroutine mxmur3(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1603
mxmu4
subroutine mxmu4(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1184
mxmfb_15
subroutine mxmfb_15(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2510
mxmur2_15
subroutine mxmur2_15(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1549
mxmur3_8
subroutine mxmur3_8(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1960
mxmf3_18
subroutine mxmf3_18(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3284
mxf12
subroutine mxf12(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:271
mxmfb_10
subroutine mxmfb_10(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2390
mxf5
subroutine mxf5(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:131
i
for i
Definition: xxt_test.m:74
mxmfb_9
subroutine mxmfb_9(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2369
mxmur3_13
subroutine mxmur3_13(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1790
mxmfb_5
subroutine mxmfb_5(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2295
mxm44
subroutine mxm44(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3515
mxmf3_2
subroutine mxmf3_2(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2940
mxf24
subroutine mxf24(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:632
mxmur2_7
subroutine mxmur2_7(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1377
mxm44_2_t
subroutine mxm44_2_t(a, m, b, k, c, n)
Definition: mxm_std.F90:3976
mxmf3_10
subroutine mxmf3_10(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3080
mxmur3_7
subroutine mxmur3_7(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1988
mxf8
subroutine mxf8(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:185
mxmf3_3
subroutine mxmf3_3(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2954
mxmur3_1
subroutine mxmur3_1(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2114
mxmf3_14
subroutine mxmf3_14(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3174
mxmur2_14
subroutine mxmur2_14(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1524
mxmfb_19
subroutine mxmfb_19(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2624
mxmf3_13
subroutine mxmf3_13(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3149
mxmf3_22
subroutine mxmf3_22(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3410
mxmf3_8
subroutine mxmf3_8(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3039
mxmfb_2
subroutine mxmfb_2(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2250
mxmfb_17
subroutine mxmfb_17(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2565
mxf22
subroutine mxf22(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:563
mxmur3_5
subroutine mxmur3_5(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2038
mxmf3_1
subroutine mxmf3_1(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2927
mxm44_2
subroutine mxm44_2(a, m, b, k, c, n)
Definition: mxm_std.F90:1098
mxmur2_10
subroutine mxmur2_10(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1434
mxmf3_9
subroutine mxmf3_9(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3059
mxmur3_3
subroutine mxmur3_3(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2080
mxf17
subroutine mxf17(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:408
mxf4
subroutine mxf4(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:115
mxf3
subroutine mxf3(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:100
mxmfb_21
subroutine mxmfb_21(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:2687
mxmur3_11
subroutine mxmur3_11(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1864
mxmf3_17
subroutine mxmf3_17(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:3255
mxmur2_5
subroutine mxmur2_5(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1344
mxmur3_15
subroutine mxmur3_15(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1708
mxf2
subroutine mxf2(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:86
mxmur3_9
subroutine mxmur3_9(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1930
mxmur3_10
subroutine mxmur3_10(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1898
mxf16
subroutine mxf16(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:373
mxf18
subroutine mxf18(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:437
mxmur2_12
subroutine mxmur2_12(a, n1, b, n2, c, n3)
Definition: mxm_std.F90:1477
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8