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
mxm88_(n,matA,mapA,m,matB,mapB,iwork,work)39 void mxm88_( 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: mxm88.
47 *
48 */
49
50 Integer me, nvecsA, nvecsB;
51 Integer i, k;
52
53 DoublePrecision **iptr;
54
55 DoublePrecision **ptrA, **ptrB;
56 DoublePrecision *scratch;
57
58 extern Integer mxmynd_();
59 extern Integer count_list();
60
61 extern void mxm88();
62
63
64
65 me = mxmynd_();
66
67 nvecsA = count_list( me, mapA, n);
68 nvecsB = count_list( me, mapB, m);
69
70 ptrA =(DoublePrecision **) work;
71 ptrB =(DoublePrecision **)( work + nvecsA + 1);
72 iptr =(DoublePrecision **)( work + nvecsA + nvecsB + 2);
73 scratch =(DoublePrecision *)( work + nvecsA + nvecsB + nvecsB + 3 );
74
75 /*
76 * Matrix A in A * B is triangular or symmetric.
77 */
78
79 k = 0;
80 nvecsA = 0;
81 for( i = 0; i < *n; i++ ) {
82 if( mapA[i] == me ) {
83 nvecsA++;
84 ptrA[nvecsA-1] = &matA[k];
85 k += *n - i;
86 }
87 }
88
89 /*
90 * Matrix B in A * B is full.
91 */
92
93 for( i = 0; i < nvecsB; i++ ) {
94 ptrB[i] = &matB[*n * i];
95 }
96
97
98 mxm88( n, ptrA, mapA, m, ptrB, mapB, iwork, scratch, iptr );
99
100 return;
101 }
102