1 /*
2 $Id$
3 *======================================================================
4 *
5 * DISCLAIMER
6 *
7 * This material was prepared as an account of work sponsored by an
8 * agency of the United States Government. Neither the United States
9 * Government nor the United States Department of Energy, nor Battelle,
10 * nor any of their employees, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
11 * ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
12 * COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT,
13 * SOFTWARE, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT
14 * INFRINGE PRIVATELY OWNED RIGHTS.
15 *
16 * ACKNOWLEDGMENT
17 *
18 * This software and its documentation were produced with Government
19 * support under Contract Number DE-AC06-76RLO-1830 awarded by the United
20 * States Department of Energy. The Government retains a paid-up
21 * non-exclusive, irrevocable worldwide license to reproduce, prepare
22 * derivative works, perform publicly and display publicly by or for the
23 * Government, including the right to distribute to other Government
24 * contractors.
25 *
26 *======================================================================
27 *
28 * -- PEIGS routine (version 2.1) --
29 * Pacific Northwest Laboratory
30 * July 28, 1995
31 *
32 *======================================================================
33 */
34 #include <stdio.h>
35 #define min(a,b) ((a) < (b) ? (a) : (b))
36
37 #include "globalp.c.h"
38
mxm5x_(n,matA,mapA,m,matB,mapB,iwork,work)39 void mxm5x_( n, matA, mapA, m, matB, mapB, iwork, work)
40 Integer *n, *mapA, *m, *mapB, *iwork;
41 DoublePrecision *matA, *matB, *work;
42 {
43 /*
44 *
45 * This is the FORTRAN wrapper for the matrix multiplication
46 * routine: mxm5x.
47 *
48 */
49
50 Integer me, nvecsA, nvecsB;
51 Integer i, k;
52
53 DoublePrecision **ptrA, **ptrB;
54 DoublePrecision *scratch;
55
56 extern Integer mxmynd_();
57 extern Integer count_list();
58
59 extern void mxm5x();
60
61
62
63 me = mxmynd_();
64
65 nvecsA = count_list( me, mapA, n);
66 nvecsB = count_list( me, mapB, m);
67
68 ptrA =(DoublePrecision **) work;
69 ptrB =(DoublePrecision **)( work + nvecsA + 1);
70
71 scratch =(DoublePrecision *)( work + nvecsA + nvecsB + 2 );
72
73 /*
74 * Matrix A in A * B is triangular or symmetric.
75 */
76
77 k = 0;
78 nvecsA = 0;
79 for( i = 0; i < *n; i++ ) {
80 if( mapA[i] == me ) {
81 nvecsA++;
82 ptrA[nvecsA-1] = &matA[k];
83 k += *n - i;
84 }
85 }
86
87 /*
88 * Matrix B in A * B is full.
89 */
90
91 for( i = 0; i < nvecsB; i++ ) {
92 ptrB[i] = &matB[*n * i];
93 }
94
95
96 mxm5x( n, ptrA, mapA, m, ptrB, mapB, iwork, scratch);
97
98 return;
99 }
100