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 /*
35 Internal PeIGS routine which computes
36
37 Q.W where Q is n1 by n2 matrix and W is an n2 by m matrix
38
39 the result is stored in Z
40
41 */
42
43
44 #include <stdio.h>
45
46 #include "globalp.c.h"
47
48 #define min(a,b) ((a) < (b) ? (a) : (b))
49
mxm25_(n1,n2,rowQ,mapQ,m,colW,mapW,colZ,iwork,work)50 void mxm25_( n1, n2, rowQ, mapQ, m, colW, mapW, colZ, iwork, work)
51 Integer *n1, *n2, *mapQ, *m, *mapW, *iwork;
52 DoublePrecision *rowQ, *colW, *colZ, *work;
53 {
54 /*
55 this is the wrapper for the matrix multiplication routine: Q.W
56
57 the matrix Q and W are assume to be in a column wrapped format
58
59 */
60
61 Integer me, nvecsQ, nvecsW, nvecsZ;
62 Integer *iscrat, i;
63
64 DoublePrecision **matQ, **matW, **matZ;
65 DoublePrecision *dscrat;
66
67 extern Integer mxmynd_();
68 extern Integer fil_mapvec_();
69 extern void mxm25();
70 extern Integer count_list();
71
72
73 me = mxmynd_();
74 nvecsQ = count_list( me, mapQ, n1);
75 nvecsW = count_list( me, mapW, m);
76 nvecsZ = nvecsW;
77
78 matQ =(DoublePrecision **) work;
79 matW =(DoublePrecision **)( work + nvecsQ);
80 matZ =(DoublePrecision **)( work + nvecsQ + nvecsW);
81
82 dscrat =(DoublePrecision *)( work + nvecsQ + nvecsW + nvecsZ );
83
84 for( i = 0; i < nvecsQ; i++ )
85 matQ[i] = &rowQ[i * *n2];
86
87 for( i = 0; i < nvecsW; i++ )
88 matW[i] = &colW[i * *n2];
89
90 for( i = 0; i < nvecsW; i++ )
91 matZ[i] = &colZ[i * *n1];
92
93 iscrat = iwork;
94
95 /*
96 use mxm*.c( depending on your data distribution )
97 */
98
99 mxm25( n1, n2, matQ, mapQ, m, matW, mapW, matZ, iscrat, dscrat );
100
101 /*
102 for( i = 0; i < nvecsW; i++ )
103 for( j = 0; j < *n1; j++ )
104 fprintf(stderr, " i = %d j = %d Z = %g \n", j, i, matZ[i][j]);
105 */
106
107 return;
108 }
109
110
111
112