1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
/* MatMul.c */
#include "../Bridge.h"
#define MYDEBUG 1
#if MYDEBUG > 0
static int count_MatMul = 0 ;
static double time_MatMul = 0.0 ;
#endif
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------------
purpose --- to compute a matrix-vector multiply y[] = C * x[]
where C is the identity, A or B (depending on *pprbtype).
*pprbtype -- problem type
*pprbtype = 1 --> vibration problem, matrix is A
*pprbtype = 2 --> buckling problem, matrix is B
*pprbtype = 3 --> matrix is identity, y[] = x[]
*pnrows -- # of rows in x[]
*pncols -- # of columns in x[]
created -- 98aug11, cca & jcp
-------------------------------------------------------------
*/
void
MatMul (
int *pnrows,
int *pncols,
double x[],
double y[],
int *pprbtype,
void *data
) {
int ncols, nent, nrows ;
#if MYDEBUG > 0
double t1, t2 ;
MARKTIME(t1) ;
count_MatMul++ ;
fprintf(stdout, "\n (%d) MatMul()", count_MatMul) ;
fflush(stdout) ;
#endif
nrows = *pnrows ;
ncols = *pncols ;
nent = nrows*ncols ;
if ( *pprbtype == 3 ) {
/*
--------------------------
... matrix is the identity
--------------------------
*/
DVcopy(nent, y, x) ;
return;
} else {
Bridge *bridge = (Bridge *) data ;
DenseMtx *X, *Y ;
double alpha[2] = {1.0, 0.0} ;
/*
---------------------------------
setup x and y as DenseMtx objects
---------------------------------
*/
X = bridge->X ;
DenseMtx_init(X, SPOOLES_REAL, 0, 0, nrows, ncols, 1, nrows);
DVcopy (nent, DenseMtx_entries(X), x);
if ( bridge->msglvl > 2 ) {
fprintf(bridge->msgFile, "\n inside MatMul, X") ;
DenseMtx_writeForHumanEye(X, bridge->msgFile) ;
fflush(bridge->msgFile) ;
}
Y = bridge->Y ;
DenseMtx_init(Y, SPOOLES_REAL, 0, 0, nrows, ncols, 1, nrows);
DenseMtx_zero(Y);
if ( *pprbtype == 1 ) {
/*
---------------------------------------
... vibration case matrix is 'm' or 'b'
---------------------------------------
*/
InpMtx_sym_mmm(bridge->B, Y, alpha, X);
} else {
/*
--------------------------------------
... buckling case matrix is 'k' or 'a'
--------------------------------------
*/
InpMtx_sym_mmm(bridge->A, Y, alpha, X);
}
if ( bridge->msglvl > 2 ) {
fprintf(bridge->msgFile, "\n inside MatMul, Y") ;
DenseMtx_writeForHumanEye(Y, bridge->msgFile) ;
fflush(bridge->msgFile) ;
}
/*
--------------------------------
copy solution into output vector
--------------------------------
*/
DVcopy (nent, y, DenseMtx_entries(Y) );
}
#if MYDEBUG > 0
MARKTIME(t2) ;
time_MatMul += t2 - t1 ;
fprintf(stdout, ", %8.3f seconds, %8.3f total time",
t2 - t1, time_MatMul) ;
fflush(stdout) ;
#endif
return ; }
/*--------------------------------------------------------------------*/
|
